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ABSTRACT 


Five  computer  programs  especially  useful  in  statistics  are  described  and  listings  ».i 
BASIC  are  given.  The  listings  are  generated  using  the  Hewlett-Packard  85  and  9845  desktop 
„jmpi!*ers.  The  programs  supply  the  values  of:  (1)  the  integral  of  the  bivariate  circular 
normal  distribution  (ND)  over  an  offset  circle;  (2)  the  integral  of  the  bivariate  elliptical  ND 
over  a  circle  centered  at  the  origin;  (3)  the  integral  of  the  bivariate  elliptical  ND  over  an 
offset  circle;  (4)  the  integral  of  the  bivariate  correlated  elliptical  ND  over  an  arbitrary 
polygon;  (5)  the  maximum  likeliho.  .  estimates,  obtained  from  quanta!  response  experi¬ 
ments,  for  the  mean  and  variance  of  a  ND. 
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I.  INTRODUCTION 


With  the  increasing  capability  of  desktops  for  scientific  computation,  the  need  arises  to  make 
available  in  BASIC  several  important  statistical  programs  which  have  been  operating  in  Fortran  on 
large  computers.  These  programs,  as  described  herein,  although  not  particularly  lengthy,  are  sophis¬ 
ticated  in  their  mathematical  and  logical  structure;  their  design  for  desktops  are  in  keeping  with  the 
high  standards  of  the  CDC  6700  mathematics  subroutine  library  maintained  at  Dahlgren,  [  14] . 

We  proceed  to  describe  in  mathematical  terms  the  statistical  functions  that  are  involved. 
Details  beyond  those  given  here  and  in  the  next  section  can  be  found  in  the  references. 

The  first  and  second  programs,  which  appear  together  in  one  package,  are  titled  CIRCV  and 
GCEF,  respectively. 

The  purpose  of  CIRCV  is  to  evaluate  P(R,  d),  where 

P(R,  D)  =  exp(-D2/2)  J  exp  (-r2/2)I0(rD)  r  dr.  (1) 

■'o 

Here  R  =  R/ax,  D  -  \/h*  +  k2/ax ,  where  R  is  the  radius  of  a  circle  C  in  the  xy-plane,  offset  a 
distance  V^h2  +  k2  from  the  origin.  l0(u)  is  the  modified  Bessel  function  of  the  first  kind  of 
zeroeth  order.  Statistically,  P  represents  the  probability  of  a  shot  falling  in  C,  under  a  bivariate 
normal  distribution  with  mean  zero  and  with  equal  standard  deviations  ox,  oy.  (5],  [6].  The 
function  P  is  called  the  Circular  Coverage  Function. 

The  objective  of  GCEF  is  to  evaluate  the  probability  function: 

F(K,  c >  =  c"  ^  exP  (" f  r2)  l»(yr2)  rdr,  (2) 

where 


0  <  c  a  oy/ox  <  1 ,  K  s.  R/ox  , 

A  s  (1  -c2)/2cJ,  B  s  (I  +  c2)/2c2  . 


(3) 


The  function  F(K,  c)  is  known  as  the  Generalized  Circular  Error  Function,  15),  [  15] . 

The  third  program,  called  ELLCV,  is  a  generalization  of  the  first  two.  It  supplies  the  value  of 
P,  where 
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Statistically,  P  represents  the  probability  of  a  shot  falling,  under  a  bivariate  normal  distribution 
with  mean  zero  and  with  standard  deviations  ax  and  ay  in  the  x  and  y  directions,  in  a  circle  C, 
centered  at  (h,  k)  with  radius  R,  i.e.. 


C:  (x-h)2  +  (y-k)2  =  R2.  (5) 

We  call  P(R,  h,  k,  ax,  ay)  the  Elliptical  Coverage  Function,  [3] ,  [4], 

ELLCV  is  slow  in  comparison  with  CIRCV  or  GCEF.  Thus,  we  also  give  listings  for  the  pro¬ 
gram  ELLCV3  which  supplies  P(R,  h,  k,  ax,  ay)  at  reduced  accuracy  but  at  roughly  one-half  the 
computing  time  pet  case  of  ELLCV. 

The  fourth  program  is  named  POLYCV,  it  makes  available  the  values  of  P(I1)  and  A(I1)  or 
P( A 1 ),  where 


,,  ,  (l-c2)-1'2 

Z(x-!')  '  2,SxSy-  exp 


and  A(n)  represents  the  area  of  IL  Here  n  denotes  an  arbitrary  polygon  in  the  xy-plane  which  is 
defined  by  the  coordinates  of  its  vertices  (xh  yj),  i  =  1 ,  2  •  •  •,  N.  A1  denotes  a  semi-infinite  angular 
region  in  the  plane  (see  Figures  6  and  7,  page  17).  The  integrand  of  (6),  given  by  (7),  represents  a 
correlated  bivariate  normal  density  function  with  mean  (Mx>  My)  and  covariance  matrix 


with  correlation  coefficient  c,  (id  <  1).  Details  of  the  analysis  for  POLYCV  are  given  in  [9], 
(101,  HI],  (12]. 

The  last  program  is  named  MLEQRE,  which  stands  for  “Maximum  Likelihood  Estimates  from 
Quantal  Response  Experiments.”  An  estimated  mean  ft  and  an  estimated  standard  deviation  a  of  a 


•NotMhat  (l)and  (2)  follow  from  (4)  by  transforming  (x,  y)  to  polar  coordinates  (r,  8)  and  setting  ox  *oy  for(l), 
and  h  =  k  =  0  for  (2). 
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normal  distribution  obtained  from  quantal  responses  are  supplied  by  MLEQRE.  It  also  gives  the 
associated  covariance  matrix  elements  and  a  plot  of  the  50  and  95%  confidence  ellipses.  The  esti¬ 
mates  o  are  taken  as  those  unique  values  of  the  independent  variables  u,  s,  respectively,  which 
maximize  the  likelihood  function 


where 


N  M 

f  =  n^>  n  q^). 

j  =  l 


i  =  1 


Z(v) 


=  -^==  exp  (~v2/2) ,  P(v)  =  J  Z(r)  dr , 


(9) 


Q(v)  =  f  Z(r)dr, 


(10) 


Sj  -  (as  -  u)/s ,  tj  =  (bj  -  u)/s . 


(ID 


The  a^  bj  are  input  stimuli  from  a  set  of  quantal  experiments.  They  are  called  “successes”  and 
“failures,”  respectively,  [7],  [8],  where  N  denotes  the  number  of  successes  and  M  the  number  of 
failures. 

In  the  next  section  we  discuss  each  program  and  how  to  use  it.  The  third  section  contains  the 
program  listings  in  BASIC  with  sampie  outputs. 


II.  MATHEMATICAL  DESCRIPTION  OF  PROGRAMS 


In  this  section,  the  five  programs,  introduced  in  the  previous  section,  namely 

CIRCV 

GCEF 

ELLCV  (ELLCV3) 

POLYCV 

=bRE 


MLE| 


are  described  in  mathematical  terms. 


All  of  the  programs,  except  POLYCV,  contain  a  subroutine  for  computing  the  complementary 
error  function,  erfc  (•),  to  a  preset  relaiive  accuracy  using  Cody’s  rational  functions,  [2], 


i.e.. 


j-i  /J-i 

erfc  (x)  =  »  -  x  £  P9(J-i)x2J/ 2  Q9(J  -  1  -  i)x2‘,  Q9(0)  s  1,  |x|  <  1/2,  (12) 

i-0 


j  l  - 1  /  j  ■  -> 

erfc  (x)  «  e-*J  £  R9(i  1  -  i)xj/  Jj 

i-0 


S' 

•  •Q 


S9(J1  -  1  - i)x',  S9(0)  =  1,  1/2  <x  <4,  (13) 


’  X.  -."  v  --C  »-• ; 


t  >  '■*  > 
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erfc  (x)  = 


2  r  ,  ,  n-i 

—  -W  +  -1  r 

*  X2  “ 


/J2-1 

V9(J2-i)x-2/  7  V9(J2 


- 1  -i)x-2il. 


W9(0)  =  1,  x  >  4. 


We  also  use  the  fact  that 


erfc(x)  -  2  —  erfc(-x).  (15) 

The  Cody  coefficients  are  stored  in  arrays  for  each  of  the  five  programs  as  noted  in  ( 1 2)— ( 1 4).  For 
example  in  MLEQRE  they  are  listed,  starting  with  P9(  1 ),  in  lines  225-260  of  the  HP-85  listing  (see 
page  91).  These  sets  of  coefficients  may  vary  from  one  program  to  another  depending  on  the 
accuracy  desired. 


The  subroutine  CIRCV  provides  the  value  of  P,  where 

P  =  2ff  ffexp  [“2’  ^x2  +y2'']  dx  dy’ 

c 

and  C  denotes  the  circle: 

C:(x-h)2  +  (y-k)2  =  R2 ,  (17) 

i.e.,  as  noted  earlier,  C  is  the^  circle  in  the  xy-plane  with  center  at  (h,  k)  and  radius  R.  The  normal¬ 
ized  offset  distance  from  (h,  k)  to  the  origin  is  denoted  by 

D  =  y/h2  +  k 2/cfx  .  (18) 

The  integration  of  (16)  is  carried  out  in  polar  coordinates  as  indicated  by  (1).  The  derivation  of 
( 1 )  from  ( 1 6)  is  given  in  ( 5 1 . 

Two  sets  of  recurrence  relations  are  used;  the  choice  depending  on  the  value  of  RD, 
(R  =  R/ox).  , 

For 

0  <  RD  <  7.0  , 

we  have 


(19) 

(20) 
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jgn  =  8n-l  ~^i  (y)  e‘R2/2»  SO  “  (1  —  e"R2/2) 

lk”  =  (t)  n  kn-i’  ko  *  e"°2/2 


Then  (20)  and  (21)  can  be  used  tc  obtain 
fn  =  n  +  1 


T,  =  gnk„  =(-T)4Tn-l  -  Kn 


Si  -  SI  +  Tn ,  S2  =  S2  +  Kr 


K0  =  exp  (-(R2  r  D2>/2 ] ,  T0  = 


¥  sxp  (— D2/2)  if  ^  <  5  X  10"4 


exp  (~D2/2)  -  K0  if  y  >  5X 


The  recurrence  relations  in  (22)  are  cycl'd  starting  with  n  =  0,  SI  =  T0,  S2  =  K0.  At  the  end  of 
each  cycle,  test 

(a)  n  >  (KD/y/2)  -  1 . 

Cycle  (22)  until  (a)  is  true,  then  test 

(b)  Tn  <  e . 

Continue  to  cycle  (22)  and  test  only  (b)  at  the  end  of  each  cycle  until  (b)  holds.  Then,  with 
6-decimal-digit  accuracy. 

f  P  =  SI 


1  3P  (23) 

The  function  3P/3R  is  available  at  virtually  no  additional  computation.  If  R  is  desired  for  a  fixed 
P  and  D,  then  3P/3R  can  be  used  to  find  R  by  the  Newton-Raphson  procedure. 

The  recurrence  relations  given  above,  based  on  [  1  ] ,  yield  a  slightly  faster  algorithm  than  those 
given  by  (22)  and  (23)  in  (5),  The  resulting  algorithm  and  test  (a),  as  described  above,  are  slightly 
unproved  over  those  in  ( 1 J  *. 

•We  wish  to  thank  Peter  Shugart  for  bringing  [  1 )  »o  our  attention. 
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For  large  RD  the  above  algorithm  is  inefficient.  Consequently,  the  following  algorithm,  as 
developed  in  [5] ,  [6] ,  is  used  if 

RD  >  7.0. 

Initially,  we  compute 


'x,-i  ' 


4  2  ./2RD  V? 


exp  [— (R  -  D)2/2] 


< 


(R  +  D)  1  r  \R~D\ 

M'  ~7T  vm  erfc  [~7T 


VJ 


SI  =  M,  ,  S2  =  Xj  . 


(24) 


Then  starting  with  n  =  0: 


n  =  n  +  1 


V  -  (¥)  4g 


4RD  y'2n'1 
M2n  +  1  =  iR2  -  D2|X2n  +  i 
^2n*l  =  (2n—  l)X2n+1 


M2n-! 


SI  =  SI  +  M2n  +  1  ,  S2  =  S2  +  X2n  +  ]  . 


(25) 


Iterate  (25)  until 


^2n  +  l  ^  e  • 


(26) 


When  (26)  is  satisfied,  continue  to  cycle 


n  =  n  +  1 


y  (2n-  l)2  1  y 

2n<1  "  ,  2n  4RD  A2n_1 


S2  =  32  +  X2n+1 


(27) 


until 


X2n*l  <  6  • 


(28) 


•The  quantity  X2n+j  as  given  above  differs  by  a  factor  of  1/4RD  from  X2n+I  as  given  in  (31),  (32),  (34)  of  (6) . 
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When  (28)  holds,  then  with  6  decimal-digit  accuracy 


fp  = 


j  11  +  sgn(R-D)  -  S2  -  sgn  (R  -  D)*S1 1 


< 

3P 

^9R 


R*S2 , 


where  sgn  (R  -  D)  =  0  if  R  =  D. 


(29) 


We  note  that  Cody’s  algorithm  for  erfc  (•)  is  only  needed  for  arguments  <4.  There  is  no  need 
to  consider  larger  arguments  because  if  | R  —  D|/\/T  >  4  then  P  =  1  or  0,  within  6  decimal  digits, 
depending  on  the  sign  of  (R  -  D),  (see  lines  205,  210  of  the  HP-85  listing).  In  case  it  is  desired  to 
achieve  greater  accuracy,  then  for  arguments  larger  than  4,  use  in  place  of  Mj  in  (24) 


(R  +  D\  \P- 

k  spl  )  IR-DI 


Of  course,  other  changes  would  also  be  needed,  such  as  changing  the  value  of  e  and  the  right  hand 
side  of  RD  >7. 


The  subroutine  GCEF  is  included  in  the  same  program  package  with  CIRCV,  because  the 
Generalized  Circular  Error  Function  F(K,  c)  can  be  obtained  using  P(R,  D).  Indeed,  we  have  from 
(16)  and  (17)  in  [5] 


P(R,  D)  -  P(D,  R)  ■  sgn  (R  -  D)F(K,  c) 

P(R,  D)  +  P(D,  R)  =  1  -  exp  [ — (R2  +  D2)/2]  I0(RD) . 


(30) 


But  from  (33)  below  R  -  D  >  0,  hence 

F(K,  c)  =  |2P(R,  D)  -  1  +  exp  [-(R2  +  D2)/2]10(RD)| .  (31) 

Note  that  the  product  in  (31)  of  the  exponential  and  Bessel  function  is  given  by  S2  in  the  computa¬ 
tion  of  P(R,  D). 

The  arguments  of  F(K,  c)  are  given  by 


f  K  =  R/ox  ,  (R  comes  from  (3)) * 
1  0  <  c  *  ay/ox  <  1 , 


(32) 


where  R  and  D  in  (31)  are  expressed  in  terms  of  K  and  c  by 


D  *  K 


c^0, 


(33) 
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and 


RD  =  K2 


1  -c2 


1C1 


It  can  also  be  shown  that 


||  =  (K/c)S2,  c  *  0. 

The  case  c  =  0  is  treated  separately  in  GCEF  with 

'F(K.O)  =  1  -  erfc(K/>/2) 

IIk  <k-0)  “V? S2-  (* /?< ‘~K’12) ■ 

We  note  if  c  >  1,  then  simply  redefine  K  and  c  so  that 

JK  =  R/oy  ,  (R  from  (3» , 

U  =  ox/oy  , 


(34) 


(35) 


(36) 


(37) 


i.e.,  interchange  ax  and  ay,  (see  (2)). 

The  program  input  variables  are  R,  D,  V.  The  output  variables  are  P,  S2. 

CIRCV 

Input:  R/ox,  D/ox,  1(=V),  (ax  >0). 

Output:  If  R/ax  and  D/ax  >  0,  then  P  =  P(R,  D),  S2  =  1/R  3P/3R.  If  R  and/or  D  <  0, 
then  P  =  -1  indicating  unacceptable  input. 

GCEF 

Input:  R/ox(=K!),  ay/ox(=c),  0(=V),  (ax,  ay  >0). 

Output:  If  R/ox  >  0  and  0  <  oy/ox  <  I .  then  P  =  F(K,  c)  and  S2  =  c/K  3F/3K(c  =£  0); 
S2  =  VirTf  3F/3K  if  c  =  0.  If  R/ax  <  0  or  |oy/ax  -  1/2|  >  1/2  then  P  =  -1 
indicating  unacceptable  input.  If  ay  >  ax  the  user  must  interchange  ax  and  ay. 


The  subroutine  ELLCV  provides  the  value  of 

c 


~  2iraxoy  //CXp 
c 


1 

2 


P 


dx  dy ; 


(38) 
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C  denotes  the  circle  given  by  ( 1 7).  It  is  shown  in  (31  that  (4)  can  be  expressed  in  the  form 


P  =  --4=  f  (exp  (-Xq/2)  +  exp  (~X[/2))  (erl'c  (y0)  -  erfc(y,))  tdt, 
V  (j 


(39) 


where 


X0  a  h  -  "(1-t2),  y0  = 


k  -  Rt  V2  -  t2 


Xi  a  h  +  R(1  - t2),  y,  = 


‘A 

k  +  Rt  V2  -  t- 


(40) 


R  =  R/ox ,  R  =  R/o  ,  h  =  h/ox ,  k  =  k/ov . 


Without  loss  of  generality  h  and  k  are  assumed  to  be  non-negative. 

The  average  computation  time  to  evaluate  P  from  (39)  is  an  order  of  magnitude  larger  than  the 

average  time  for  CIRCV  or  GCEF.  Hence,  a  number  of  tests  are  used  to  determine  if  P  <  e  or 

P  >  1  -  e  in  which  case  P  is  set  to  zero  or  om  ,  respectively.  Let  H2  =  h2  +  k2,  o  =  max  (ox,  ay). 

Test  #/:  If  R2  <  2eoxoy  then  P  <  e. 

Test  #2.  If  R  -  h  +  A ax  <  0  or  R  -  k  +  Aoy  <  0  then  P  <  e. 

Test  #2.  If  R  >  \/h2  +  k2  +  Al*o,  then  P  >  1  -  e. 

Test  #4:  If  H2  >  R2  and  if 

R2  exp  j— i  ((H-R)/o)2J  <  2sotOy,  then  P  <  e. 


The  value  of  A!  is  chosen  so  that  E  (see  Figure  1)  contains  1  -e  of  the  distribution;  A  is  chosen  in 
a  similar  way  with  E  replaced  by  a  rectangle  centered  at  the  origin  with  sides  of  length  2A ox  and 
2Aoy  along  the  x  and  y  axes,  respectively,  (See  (50)  and  (5 1 )). 

Test  #1  follows  by  taking  H  =  0  and  considering  small  R.  Tests  #2  and  #3  are  covered  in  (3) . 
Test  #4  follows  by  using  the  fact  that 


exp  <-■ 


H 


x 


<  exp 


2o2 


(x2  +'y2) 


Tests  #5  and  #6,  given  below,  are  more  subtle.  We  discuss  test  #5  first. 
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Hence,  if 


Consider  Figure  ( 1 ).  A  tangent  line  L  is  drawn 
at  point  B  where  the  ray  M,  from  the  origin  to  (h,  k), 
intersects  C;  consider  also  a  line  L  with  the  same  slope 
as  L  tangent  to  E,  the  ellipse  with  center  at  the  origin 
with  major  and  minor  axes  2Aioy,  2Alox.  The 
distance  from  the  origin  to  B  is  given  by 

d(L)  =  H  -  R  >  0  (41) 

and  the  normal  distance  from  L  to  the  origin  is  given 
by 

d<L)  =  Alox[(h2  +  a2k2)/H2l1/2 ,  o  s  0y/ox . 

(42) 


d(L)  >  d(L) 


(43) 


then  P  <  e.  This  test  can  be  formulated  without  using  square  roots  and  is  called  Test  #S. 


Test  #5:  If 


y  =  H2  -  R2  -  Al2ox2 


>  0 


(44) 


and 

y2  >  4R2A1  (45) 


then  I’  <  e. 

For  test  #6  consider  Figure  (2).  Here  L  denotes  the  tar.jent  line  to  E  at  D  where  ray  M 
intersects  E  at  D;  L  denotes  the  line  tangent  to  C  which  intersects  M  and  which  is  parallel  to  L. 
In  this  case 

T  s  d(L)  -  d(L)  *  H{l  -  R  v^/D  -  Al/y^}  (46) 

where 

F  =  o2h2  +  k2 ,  D  =  h2  +  k*.  R  s  R/oy. 

If  T  >  0,  then  P  <  e.  Without  using  square  roots,  we  have: 
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Test  #6:  If 


y  =  D  -  R2F/D  -  Al2  >  0, 


and  if 


then  P  <  e. 


y2  >  4A12R2F/D, 


If  none  of  the  above  tests  are  applicable,  then  a 
Gaussian  numerical  integration  is  used  to  evaluate  P 
from  (39).  In  this  case  it  is  advantageous,  if  possible, 
to  reduce  the  interval  of  integration.  By  methods  similar 
to  those  described  on  pp  7-9  of  [3],  an  interval  of 
integration  [e0,  ej  C  [0,  I]  is  determined  such  that 


Figure  2.  Test  #6 


(e0,  ejl  =  min  {[e0,  e,],  [e0,  e,]}  , 


where 


,  /  R  —  h  ~  A3  . .  .  _ _ 

Y - g -  if  h  +  A3  <  R, 

0  if  h  +  A3>R. 

A3-  if  0  <  h  ~  A3  <  R , 

1  if  h<  A3, 

0  if  h  -  A3  >  R. 


R  -  k  -  A3 


if  k  +  A3  <  R, 


0  if  k  +  A3  >  R . 


R_  k  +  A3  if  0<k_A3<R> 

K 

1  if  k  <  A3. 

0  if  k-A3>R. 


If  [e,  -  e0]  >  (ei  -  e0] ,  then  ax  is  interchanged  with  ay  and  h  with  k.  This  is  equivalent  to 
reversing  the  order  of  integration  in  (39). 

In  order  to  retain  efficiency  in  the  Gaussian  quadrature  evaluation  of  (39),  an  empirical  func¬ 
tion  N  is  used  to  specify  apriori,  the  order  of  the  Gaussian  process,  0(G),  to  use. ,  It  is  given  by 


-e, 
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.s 


I 


I 


I 

B 


. 

■% 

I 


For  example  if,  from  (52),  N  =  1 .3 1,  then  0(G)  =16.  This  means  16  Gaussian  abscissae  and  weights 
are  used  to  evaluate  the  right  hand  side  of  (39).  The  abscissae,  starting  with  0(G)  =  6  are  stored  in 
array  X(i),  and  the  corresponding  weights  are  stored  in  array  Y(i).  Since  the  abscissae  and  weights 
have  certain  symmetry  properties  about  zero  only  half  of  them  are  actually  stored.  Thus  the  right 
hand  side  of  (39)  is  approximated  to  within  2e  by 


1 

i 

„  1  (e\  -  eo\ 

M 

V 

P-y /SFV  2  ) 

0(G).  w0  2  0, 


(54) 


i  *  -  M/2 


where 


f  Wj  =  /th  Gaussian  weight,  w;  =  wH 
(ei  ~  e0) 


(55) 


= 


(1  +Xj)  +  e0,  Xj  =  /th  Gaussian  abscissa ,  xt  =  -x_j 


(56) 


Pi  2  {erfc  [y0(tj)l  -  erfc  [yi^tj)]}  , 

\  fi  =  {exp  [-Xg(tj)/2]  +  exp  [-Xf(ti)/:]}  ,  (See  (40)). 

Further  reductions  in  computing  time  can  often  be  realized  by  the  following: 

(a)  In  the  determination  of  [e,  -  c0) ,  it  can  be  shown  that  if  h  >  A  and  if  R  >  h  -  A  Or 
R  >  h  Jr  A  then  the  second  exponential  in  (39)  is  negligible  and  can  be  dropped.  Thus 
a  variable  H5  is  introduced  and  set  to  one  if  the  exponential  has  been  dropped. 


(b)  An  exponential  in  (39)  can  also  be  dropped  if  the  absolute  value  of  its  argument  exceeds 

>r  ELLCV) 

>r  ELLCV3).  1  * 

This  feature  is  not  included  in  ELLCV  or  ELLCV3  at  present. 


28  “  ,og [ *  2  °)  V-}+  I™ 


(c)  In  case  k  ~  0,  the  quantity  Pj  in  (56)  is  replaced  by 

p,  =  2{l  -  erfc  [yj(tj)l} 


(58) 


Thus,  only  one  instead  of  two  erfc  functions  are  required  for  each  i.  (If  ax  -  oy,  then  by 
circular  symmetry,  h  can  be  replaced  by  v  h2  +  k2  and  k  by  zero.) 

(d)  The  argument  Yo(y r )  of  erfc  (•)  is  a  decreasing  (increasing)  function  of  t.  Hence  if 

•  £ -A2(y,(T)  >  A2)  then  P(  a  2(1  -  ei)  (erfc  (yj(Dl  a  C])  for  all  t  >T.  Two 

meters  are  introduced  to  take  advantage  of  this  situation  when  it  occurs.  Variable  Z 
is  set  to  one  for  the  smallest  i  (say  =j)  on  [-M/2,  M/21  for  which  y0(tj)  <-A2  so  that 
Pi  is  replaced  by  2(1  -  e,)  for  all  ti  >  tj.  Similarly  a  variable  Z3  is  set  to  One  for  the 
smallest  i  (say  =j)  on  (-M/2,  M/2  J  for  which  y,  (tj)>A2  so  that  erfc  [y,  (tj)j  is  replaced 
by  e,  for  all  tj  >  ti.  Note: ,  t,  <  ti+.i  by  the  way  the  Gaussian  abscissae  are  ordered. 
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In  addition  to  ELLCV  which  gives  P  to  at  least  6-decimal-digit  accuracy,  a  listing  is  also  included 
for  the  program  ELLCV3  which  gives  P  to  at  least  3-decimal-digits.  The  latter  differs  from  ELLCV 
in  the  alignment  of  0(G)  by  (53),  by  using  fewer  Cody  coefficients  to  corr.pi.te  trfc  (•)  to  less 
accuracy,  in  the  values  of  the  constants  A,  Aj,  A2,  A3,  £j. 

The  values  of  these  parameters  are  given  in  the  table  below. 


ELLCV,  e  =  5 (-7) 

ELLCV3,  e  =  5 (-4) 

A 

4.892 

3.291 

A1 

5.387 

3.89895 

A2 

3.8775 

2.898 

A3 

5.16 

3.70 

1.04C-8) 

1.5(-5) 

The  values  of  O(G)  forELLCV3  are  determined  by  N  from  (52)  and  the  following  inequalities: 

(  N  >  2  -  0(G)  =  8 

I  0.675  <  N  <  2  -  O(G)  =  6 
I  0.5  <  N  <  0.675  0(G)  =  4 

l  N  <  0.5  ■*  O(G)  *  3. 


_The  inputs  to  these  subroutines  are:  R,h,  k,  ex,  ay.  The  output  is  designated  by  P.  The  value 
of  R  must  be  non-negative;  ox  and  oy  must  be  positive.  Tests  are  not  included  for  these  require¬ 
ments;  no  error  parameter  is  used. 

Accuracy:  P  is  given  approximately  to  within  2e  -  10~6  (or  10* 3  for  ELLCV3). 

Note:  Constraints  imposed  on  h,  k,  ajay  in  [3]  are  no  longer  necessary. 

The  subroutine  POLYCV  supplies  the  value  of  P(IT)  as  expressed  by  the  double  integral  in  (6). 
The  evaluation  of  (6)  is  simplified  by  using  the  transformation 


Then 

p(n)  *  ~  JfcxP  [~2  (*2  +  y2)J  d*dy  (6i) 

n 

where  II  of  (6)  is  transformed  to  II  of  (6J)  by  (60).  We  call  lithe  “transformed  polygon.”  The 
program  determines  the  vertices  (xs,  yj)  of  II  from  (60)  and  evaluates  P(n)  from  (61).  POLYCV  can 
also  be  used  to  evaluate  the  double  integral  in  (6)  over  a  semi-infinite  angular  region.  This  will  be 
discussed  further  below. 
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It  is  important  in  order  to  use  POLYCV  to  understand  how  n  must  be  specified.  We  say  n  is 
positively  oriented  (PO)  if  it  is  a  simple  polygon  or  the  limit  of  a  sequence  of  simple  polygons  as 
defined  on  page  9  of  [11],  and  if  its  vertices  are  ordered  so  that  the  area  of  n  is  on  one’s  left  as  the 
segments  of  n  are  traversed  continuously  in  the  order  the  vertices  are  given.  If  the  area  is  on  the 
right,  n  is  said  to  be  negatively  oriented  (NO)  and  P  will  be  negative.  In  case  II  has  vertices  which 
occur  more  than  once  both  (PO)  and  (NO)  regions  can  occur.  Self-intersecting  (SI)  polygons,  as 
described  on  pages  13-17  [111,  can  also  be  handled  However  F  and/or  the  area  can  be  negative. 
The  interpretation  of  such  results  is  left  to  the  user. 

Two  examples,  to  help  clarify  these  ideas,  are  given  in  Figures  3  and  4  below.  In  Figure  3 
we  have  an  example  of  a  simple  polygon  which  is  PO.  The  probability  P  is  found  over  the  cross- 
hatched  region.  In  Figure  4,  a  polygon  is  shown  with  PO  and  NO  regions.  The  probability  is  found 
again  over  the  cross-hatched  areas. 


The  description  given  above  for  specifying  n  is  adequate  for  most  applications.  When  prescribing  a 
completely  arbitrary  polygon,  all  the  vertices,  points  where  two  segments  cross,  and  initial  and 
terminal  points  of  overlapping  segments  must  be  numbered.  In  certain  situations  some  of  these 
points  may  not  be  necessary  as  shown  by  the  example  on  page  24  of  [111.  However,  when  in 
doubt,  II  should  be  numbered  as  just  described.  More  details  and  examples  are  given  on  page  14 
of [  1 11. 

It  is  shown  in  [  1 1  ] ,  [  1 21  that  P  over  n  is  given  by, 

N 

POT)  »  W  -  £  P(AS),  (62) 

i-  1 

where  N  is  the  number  of  points  specifying  II,  and 


'  W  =  £2/2jt,  (See  (75)  for  12), 

P(Aj)  =  r—  j'f  exp  (V  +y2)J  dxdy. 

'•  A( 


(63) 
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Here  A;  denotes  a  semi-infinite  exterior  angular  region  of  fl.  In  Figure  5,  the  four  exterior  angular 
regions  are  shown  for  a  simple  polygon  with  N  =  4.  For  i®  1, ....  N,  A4  is  formed  by  extending  the 
side  (i  -  1,  i)  from  (i)  to  00  in  the  direction  from  (i  -  1)  to  (i).  Note  (C)  =  (N).  Similaily  side 
(i,  i  +  1)  is  extended  from  (i  +  1)  to  00  in  the  direction  (i)  to  (i  +  1),  where  (N  +  1)  =  (1).  The 
angle  of  Aj,  with  its  vertex  at  (i),  is  measured  as  positive  in  the  counterclockwise  direction  and  as 
negative  in  the  clockwise  direction. 


In  Figure  5,  angular  regions  At,  A2,  A4  have  positive  measures  with  P(A,),  P(A2),  P(A4)  positive, 
angular  region  A3  has  negative  measure  with  P(A3)  <  0. 


The  sum  of  the  angular  measures  of  the  Ajt  i  =  1, 2,  ....  N  appears  in  (63)  and  is  denoted  by  !2. 


In  order  to  evaluate  P(Aj),  advantage  is  taken  of  the  circular  symmetry  of  the  integrand  in 
(63).  A  semi-infinite  straight  line  L  is  introduced  extending  from  the  origin  to  Vj,  the  vertex  of  Ai( 
to  00 .  Then  transforming  the  integration  variables  in  (63.)  to  polar  coordinates  (r,  0)  at  Vj.  with 
6  >  0  when  measured  counterclockwise  from  L  about  Vj,  we  have  as  derived  on  pp  2-3  of  (9). 


P(Aj) 


j  (R2  +  2rR  cos  6  + 


r  dr  dd , 


(64) 


where  R  is  the  distance  from  the  origin  to  v,,  and  02  -  Q\  *  A0  is  the  angular  measure  of  At.  Note  if 
Vj  is  at  (0,  0),  so  that  R  *  0,  then  PfAj)  =  A0/2ir. 

Using  the  fact  that 


f  e"r2/2e“,R  co,0  r  dr  =  1  -  2u  erfc  (u)/*(u),  (65) 

4 
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where 


u  =  — ■=  cos  9 ,  z(u) 

V  2 


^  -OO 

=  — p=exp  (-U2),  erfc(u)  =  f  z(t)  dt, 

V*7  -l 


we  have 


'  A9  i  r 3* 

p'*l,  u 


[erfc  (u)/z(u)]  dd  ,  -ir<A9<  n. 


P(Aj)  =  e'R2/2 


The  evaluation  of  the  right  hand  side  of  (67)  is  achieved  by  using  a  minimax  polynomial  approxima¬ 
tion  on  [0,  C(S)],  i.e.,  given  a  5  >  0,  a  set  of  real  numbers  {Ulk}  are  found  for  a  least  positive 
integer  K  such  that 

K-l  » 

erfc (u)  -  z(u)  £  0 1 K _ kuk  <  -~=S,  0  <  u  <  C( 6).  (68) 

k  =  o  * w 

Given  a  value  of  6  >  0,  the  constant  C(6)  is  determined  so  that  the  value  of  P(A)  in  (63)  is  negli¬ 
gible  when  A  is  the  angular  region  with  C(6)  =  R.  with  vertex  on  the  positive  x-axis  and  with 
92  =  ~9i  =W2.  For5  =  5  X  lO'10,  C(5)  =  6.2,  and  K=  15. 

The  coefficients  {Ulk}  for  5  =  5  X  10"10  are  given  in  the  HP  9845  POLYCV  listing  in  lines 
205-275.  Coefficients  for  some  other  5’s  are  given  on  page  E-4  of  [  1 1  ] . 

With  (68)  it  is  not  difficult  to  show  that,  within  S/s/ir, 


e-R2/2  A9i 

P(At)  =  v  V 


'a9{  W 

2  Zj  UIK-kJk  +  t  » 

k-0 


where 


Ik  =  (^j  cos*9d9,  |0,|  <  j, 

J  Jk*l  ■  7T“f  {[h2g2  “  Ml]  + 


|02l  <f. 


*  ¥  '  J 


Jq  =  A5 ,  Jj  =  hj  —  hj 


gj  =-~cos0j,  hj  =-^=s:n0j,  j  =  1,  2, 


R2  =  x2  +  y2  (vertex  of  Aj  at  (x,  y)). 


(73) 
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Since  u  >  0  in  (68),  this  requires  the  constraints  on  6X,  02  given  in  (70).  For  rr/2  <  0  <  n,  we  use, 
in  addition  to  (69), 


P(A[R,0])  =  j  erfc 


P(A[R,  7T  —  0])  , 


(74) 


where  A[R,  0]  denotes  an  angular  region  with  its  vertex  a  distance  R  from  the  origin  and  with 
angular  measure  d,  where  one  side  of  A[R,  9]  is  formed  by  the  line  L  described  above.  See 
pp  13-14  of  [9]  for  more  details. 

We  note  that  since 


N 

n  =  £  Aflj, 

i*l 


(75) 


and  A0j  already  occurs  in  (69),  lPtic  additional  computing  is  necessary  to  obtain  W  in  (62). 

Since  most  of  the  programming  is  already  available,  POLYCV  is  also  designed  to  yield  the 
normal  probability  P(A1)  over  a  single  semi-infinite  angular  region  Al.*  However  for  polygons  each 
of  the  angular  measures  A0;  in  (75)  is  in  the  interval  (-ir,  it] ,  whereas  for  the  single  angular  region 
case  the  angular  measure  of  Al  is  always  in  the  interval  [0,  2i r).**  The  angular  region  A 1  is  specified 
by  three  points.  The  first  point  is  always  at  the  vertex  of  A 1  the  second  and  third  points  are  taken 
such  that  A0  of  A 1  is  >0.  Figures  6  and  7  show  typical  angular  regions. 

In  addition  to  P(=P(n)),  POLCV  also  supplies  as  output  the  area  of  II,  A(fl),  where 


AOI)  =  SxSy[l  -c2]1/2A(II), 

N 

A(n)  =  H  xi<yi*i  yo  =  ys>  yn+i  =  yi- 


i-  1 


*P(A  1 )  corresponds  to  a  Bivariate  Normal  Integral.  See  [9,  p.  18], 

••Since  A0,  in  computing  P(A1),  is  always  >0,  P(A1)  is  always  non-negative. 


(76) 

(77) 
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In  (77),  the  x,,  y;  refer  to  the  coordinates  of  the  /th  vertex  of  the  transformed  poi/gon  n  as 
obtained  from  (60).  For  PO  (NO)  specified  polygons  A(II)  is  always  >(<)0.  Thus,  as  indicated 
earlier  for  arbitrary  polygons,  one  should  call  A(H)  the  “signed”  area  of  II. 

Tne  input  for  POLYCV  is  specified  in  data  statements  accordingly: 


HP-9845:  P8,  P9,  Mx,  My,  C,  Sx,  Sy,  N,  Xj,  y3,  x2,  y2,  x*,  yK, 

HP-85:  P8,  P9,  Ml,  M2,  C,  SI,  S2,  N,  x3,  y^  x2,  y2,  ....  xK,  yK. 


For  the  HP-85,  Ml  =  Mx,  M2  -  My,  SI  =  Sx:  S2  =  Sy. 

If  P8  =  0,  then  the  vertex  coordinates  Xj,  y;,  i  =  I  to  K  are  stored  in  data  statements  im¬ 
mediately  following  the  initial  data  statement  above.  POLYCV  stores  them  in  arrays  X.(*),  Y(*). 
If  P8  #  0,  it  is  assumed  the  vertex  coordinates  are  already  stored  in  arrays  X(*),  Y(*)  which  is 
convenient  if  the  vertices  are  macnine  generated. 


P9:  Determines  the  output  desired., 

P9  =  0.  No  listing  of  the  xjt  y4 ;  no  piot  of  n. 

P9  S*  1 .  A  listing  of  the  Xj,  ys  is  printed. 

P9  >  1  or  P9  <  0.  A  plot  cf  n  or  angular  region  A1  in  the  xy-plane,  depending  on  N,  is 
given  on  the  CRT  and  dumped  onto  the  printer. 

N:  If  N  >  3,  then  K  is  set  tr  N  by  POLYCV  and  P(II)  is  found. 

If  N  -  1,  then  K  is  set  to  3  and  P(A1)  is  found. 

Xj,  yj,  i  =  1,  2,  ....  K:  (Xj,  y j)  denotes  the  /th  point  specifying  a  polygon  n  or  an  angular  region 

Al.  If  n  is  simple  the  points  should  be  ordered  so  that  n  is  positively 

oriented.  In  the  case  of  an  angular  region  AI(N=  1),  (Xj,  yj)  locates  the 

vertex  of  Al,  points  (x2,  y2)  and  (x2,  y 3)  are  ordered  so  that  one  rotates 

from  (x2,  y2)  to  (x3,  y3)  in  a  counterclockwise  direction  about  (Xj,  yt). 

Four  output  quantities  are  always  pnnted.  They  are  P,  A,  W,  1 1. 

P:  Contains  the  value  of  P(II)  if  N  >  3  or  P(A1)  if  N  =  1. 

A:  Contains  the  area  of  II  if  N  >  3  or  is  set  to  0  if  N  =  1. 

W:  Denotes  the  “winding  number  of  II,”  see  page  18  of  (111.  For  a  simple  (PO)  polygon  it  is 

always  one.  It  is  set  to  zero  if  N  3  1. 

II:  If  II  =  0  or  2  output  is  acceptable,  11  =  2  indicates  that  two  or  more  consecutive  sides  of  n 
overlap.  If  II  =  1,  then  angular  region  Al  with  N  =  1  is  not  well-defined,  i.e.,  the  vertex  of 

Al  and  at  least  one  of  the  other  two  points  specifying  Al  are  too  close  to  each  other.  If 

II  =*  —  1  then  Al  may  not  be  well-defined -the  angular  measure  of  Al  is  close  to  Oor  2 n.  A 
value  for  P  is  given.  If  II  =3,  then  c,  the  correlation  coefficient,  does  not  satisfy  c2  <  1  and 
is  unacceptable. 
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The  routine  is  presently  set  to  yield  P(II)  or  P(A1)  to  approximately  9-decimal-digit  accuracy.  The 
computing  time  can  be  reduced  significantly  by  requiring  less  accuracy.  It  is  not  difficult  to  modify 
the  program  to  do  this.  The  necessary  changes  are  indicated  in  [  1 1  ] . 

The  final  program  vfLEQRE  provides,  from  quantal  response  experiments,  the  maximum  likeli¬ 
hood  estimates  n,  a  for  the  mean  ii0  and  the  standard  deviation  'Tq  of  a  normal  distribution.  It  also 
makes  available  the  covariance  matrix  elements  and  a  plot  of  elliptical  confidence  regions. 

MLEQRE  is  based  on  the  development  given  in  [7],  [8].  It  uses  independent  variables  a  and 
P  instead  of  u  and  s  of  (9K 1 1 ),  where 

a  =  u/s,  p  =  1/s  >  0.  (78) 

It  achieves  the  maximization  of  F,  where  F  is  given  by  (9),  by  maximizing  the  logarithm  of  F  over  a 
and  p,  where 

N  M 

L (a,P)  s  £nF  =  £  q  to  PCs,)  +  £  djfnQftj),  (79) 

i* 1  j=l 

with  Z(  ),  P(  ),  Q(  )  defined  in  ( 10)  and  with 

sj  *  /3a*  -  a,  tj  =  0bj  -  a 
c,  =  the  number  of  times  a*  occurs. 

dj  =  the  number  of  times  bj  occurs.  (80) 

"  - 1  ci*  m  =  £  dj 

i»  1  j»  1 

The  aj  and  bj  are  input  stimuli  associated  with  successful  and  unsuccessful  tests,  respectively. 
Consequently,  the  a(  are  called  “successes”  ahd  the  bj  are  called  “failures.”  They  may  take  any 
real  values,  with  R  denoting  the  total  number  of  successes  and  M  the  total  number  of  failures.  In 
order  to  take  advantage  of  the  situation  where  repeated  values  of  the  aj  and/or  bj  occur,  (79)  is 
written  in  terms  of  N  and  M,  rather  than  R  and  M,  where  N  and  M  denote  the  number,  of  different 
aj  and  bj,  respectively.  This  is  an  important  feature  if  repeated  values  occur,  because,  by  (9)  and 
(79),_fhe  computation  time  to  determine  n  and  a  will  be  proportional  to  N  +  M  rather  than 
fl  +  M. 

The  values  of  a  and  P  that  maximize  L  (and  also  F)  ate  denoted  by  A  and  B,  respectively. 
Hence 

A  =  ii/o,  B  =  1  }o>  0,  (81) 
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The  maximization  is  achieved  by  using  the  Newton-Raphson  (N-R)  procedure  in  two  independent 
variables  Initial  estimates  for  A  and  B,  which  are  required  for  (N-R),  are  denoted  by  AO  and  BO. 
Either  they  are  supplied  by  the  user,  or  with 

2-1  =  W  C‘ai  ’  ^2  =  1  £  djbj  ’  £3  =  (N  +  m/^J  °iar  +  £  djb^ 

1=1  j=i  \ 1 = 1  j=i  / 


MLEQRE  uses 


TbO=  J Y 1  -  r  —  1  -  (NS,  +M12) 

1  1^-3  Un  +  m)v  1  v 


:)lf  >  0 


(82) 


AO  =  j  (ll  +  S2)  •  BO. 


The  (N-R)  corrections  Dl,  D2*,  beginning  with  corrections  to  AO,  BO,  are  found  by  solving  two 
linear  equations  (see  (55)  of  [7]  and/or  (4.1)  of  [8]  with  their  coefficients  expressed  in  terms  of 
the  first  and  second  partial  derivatives  of  L  with  respect  to  a  and  0  (see  (1 1 5)— ( 1 19)  of  [7]  and/or 
(3.6M3. 10)  of  [8] ). "For  completeness,  we  give  these  relationships  here  again. 


"(Dl)Laa  +  (D2)La0  =  -L* 


(N-R)  Equations,  Laa  = 


,(D1)LU<J  +  (D2)Lp0  =  -Le 
f  Dl  =  a.j>l^-LaL^)lA,  D2  =  (LQLa^-L^Laa)/A 


= 

da2 


X 


A  —  LqqLpp  “  L q,^  0 


M 


La  “  djYj /Qi  “X  CiXj/Pj,  yj  =  Z(tj),  Xj  s  Z(Sj) 


(83) 


(84) 


(85) 


1*  1 


N 

h  =  Z  ciWpi  ~  £  diWQi 

i =  1  j*l 


(86) 


M  N  ' 

Laa  =  -£  dj(yj/Qj)[(yj/Qj)  -  tj]  -  £  +  s,]  <0  (87) 

j-i  ^ 


i«  l 


*Some  quantities  such  as  Dl  and  D2  are  used  both  as  real  variables  and  also  as  BASIC  variables  which  contain  the 
corresponding  real  variable  values.  It  should  be  clear  from  the  context  which  is  intended. 
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N  M 

Ltt0  =  £  ciai(xi/Pi)t(xi/Pi)  +  Sj]  +2]  djbj(yj/Qj)[(yj/Qj)  -  tj]  (88) 

i=l  j=l 

M  N 

Lfifi  =  -£  djbfCyj/QjMCyj/Qj)  -  tj]  -£  qafai/PjHUj/Pj)  +  s,]  <  0.  (89) 

j=i  ,  i=i 

The  equation  for  a  confidence  ellipse  in  the  us-plane  with  center  at  (/i,  u)  is  given  by 

Auu(u-/a)2  +  2Au,(u-^)(s~o)  +  A„(s-o)2  =  x?-7>  (90) 

where  xj-y  is  obtained  from  a  chi-squared  table  with  two  degrees  of  freedom;  y  denotes  the  prob¬ 
ability  that  the  ellipse  contains  (/i0,  O0).  For  y  =  .5,  x2s  =  1-39,  and  for  y  ®  .95,  x2os  =  5.99.  See 
page  42  of  [7]  for  other  values.  The  coefficients  in  (90),  which  make  up  the  elements  of  the  inverse 
of  the  covariance  matrix,  are  available  directly  from 


o2Auu  =  £  CifXi/PiKXi/Qi)  +  £  djty/PjXyj/Qj) 

i=l  j“l 

N  M 

°2K,  =  2]  CjSiCxi/PiKxj/Qi)  +  £  djtj(yj/Pj)(yj/Qj) 


i=  1 
N 

E 

i“  1 

N 

E 

i=  1 


(91) 


j  =  i 

M 


a2AM  =  21  Cj^Cxj/Pt) (xj/ Qt )  +  2]  djt2(yi/Pj)(yi/Qj), 
*“  i‘i 


which  are  evaluated  at  a  =  n/a  and  $=  l/o;  it  is  recalled  from  (78),  (80)  that 

'sj  ■  aj/J  -  a 

tj  =  bj/3  -  a,  as  u/s,  /3  =  ,1/s  >  0. 

Derivations  of  (90)  and  (91)  are  given  in  [7] . 

The  evaluation  of  the  quantities  in  (84M90)  requires  an  efficient  and  high  precision  subroutine 
to  compute 


Y  3  y/Q,  Y1  H  (y/Q)[y/Q  - 1]  ,  Y2  m  (y/Q)(y/P), 


(92) 


where  the  subscript  j  has  been  dropped.  It  is  easy  to  show  that  the  corresponding  quantities  in 
terms  of  (Xj/Pj)  can  be  obtained  from  the  same  subroutine  by  changing  the  sign  of  the  input  argu¬ 
ment,  i.e.,  changing  Sj  to  -Sj. 
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Cody’s  rational  approximations  for  the  complementary  error  function  as  given  by  ( 1 2>—<  15) 
are  fundamental  to  the  evaluation  of  the  quantities  in  (92).  Using  J  =  4,  J1  =6  and  J2  =  4  in 
( 1 2)-(  1 5)  yields  a  minimum  of  eleven  and  one-half  significant  digits  of  accuracy  for  the  comple¬ 
mentary  error  function,  erfc(-)- 


For  completeness,  we  give  the  expressions  used  for  Y,  Yl,  Y2,  where  we  use  the  facts  that 

y(t)  =  x(t)  =  x(  t)  =  y(-t) 

(93) 

P(t)  =  Q(-t). 

Let  t  =  tj  or  (— Sj),  K1  =  t/y/T,  C2  =  y/2/ir,  E  =  exp  (-K1 2).  Also  let 

2  =  2n/Zd, 


where  LN  denotes  the  numerator  sum  in  (12),  (13)  or  (14)  and  2D  denotes  the  denominator  sum. 
For  example  if  1/2  <  K1  <  4,  then,  referring  to  (13), 


R9(J1  -‘XKIV- 


Now  for: 

|K1|  <  1/2 

r  Y  =  (C2)E/erfc(Kl) 

Yl  =  Y(Y-t) 

_  Y2  =  (C2)EY/(1  +K1  2) 
-4  <  K1  <  —1/2 

rY  =  (C2)F./(2  -  E  2) 

Yl  =  Y(Y-t) 
lY2  =  (C2)Y/2 
1/2  <  K1  <  4 
f  Y  =  C2/2 
Yl  =  Y(Y-t) 

Y2  »  (C2)EY/(2-E2) 


(94) 


(95) 


(96) 
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-5.5  <  K1  <  -4 

f  Y  =  (C2)E/[2  -erfc(Kl)]  ,  K1  =  |K1| 


<  Y1  =  Y(Y-t) 


LY2  =(C2)K1  Y/(1/v/¥+2/K12) 

K1  >  4 

"Y  «  (C2)Kl/[(l/>/»)  +  S/Kl2] 

<  Yl  =  -(2/v/i)2/[(l/^)  + 2/K12]2 
LY2  »  (C2)EY/[2  -ferfc(Kl)] 


(97) 


(98) 


K1  <  -5.5 

Y  =  Yl  =  Y2  =  0  (Y  <  3  X  10"14)  (99) 

The  subroutine  for  computing  Y,  Yl,  Y2  using  the  above  begins  on  line  1215  for  the  HP  9845 
program  and  on  line  695  for  the  HP-85  program.  It  is  the  core  of  MLEQRE. 

There  are  numerous  BASIC  variables  in  MLEQRE  which  are  pertinent;  we  list  them  here  for 
convenience. 


A(i) :  cpntains  the  /th  listed  success  value  ai,  (1  <  i  <  N). 

B(j) :  contains  the  /th  listed  failure  value  bj  (1  <  j  <  M). 

C(i) :  contains  cj5  the  number  of  aj  at  a  fixed  i. 

EXJ) :  contains  dj,  the  number  of  bj  at  a  fixed  j. 

N :  contains  the  number  of  different  aj. 

M  .  contains  the  number  of  different  bj. 

AO :  contains  the  current  o(*  u/s)  value. 

BO:  contains  the  current  0(=  1/s)  value. 

U :  contains  \i. 

S:  contains  a. 

A 1 :  contains  the  initial  estimate  for  a  upon  EXIT 

B 1 :  contains  the  initial  estimate  for  0  upo  i  EXIT 

f  If  LO  contains  one,  then  c(  and  d(  are  one  for  all  i  and  j. 
LO:  < 


.If  LO  does  not  contain  one,  then  some  c(  and/or  dj  may  be  larger  than  one. 
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D1 :  contains  the  current  (N-R)  correction  D1  to  AO. 

D2  :  contains  the  current  (N-R)  correction  D2  to  BO. 

LI :  contains  La(=  3L/3a) 

L2 :  contains 
L3 :  contains  Laa 
L4 :  contains  La^ 

L5 :  contains  Lgp 

L6,  L7,  L8 :  contain  the  confidence  ellipse  coefficients  Auu,  AUJ,  An,  respectively 

L2,  L3,  L4 :  change  their  contents,  after  completion  of  the  (N-R)  procedure,  to  the  covariance 
matrix  elements  Auu,  Au‘,  A,s,  respectively. 

Y :  contains  either  (x-J Pj)  or  (yj/Qj) 

Y1 :  contains  either  (xi/Pi)[(xi/Pi)  +  Sj]  or  (yj/Qj) [(yj/Qj)  -  tjl 

Y2 :  contains  either  (xi/Pi)(xi/Qi)  or  (yj/Pj)(yj/Qj) 

Z :  contains  0.  1,2  and  is  used  to  signal  that  one  cycle  of  (N-R)  remains  to  be  carried  out. 
This  allows  Y2  to  be  computed,  for  use  in  the  evaluation  of  the  confidence  ellipse 
coefficients  in  (90),  only  on  the  last  (N-R)  iteration. 


The  input  data  for  MLEQRE  is  stored  in  data  statements.  If  LO  #  1 ,  then  data  is  stored 
sequentially  in  the  following  variables:  N,  M,  LO,  AO,  BO,  P8,  C(l),  A(l),  C(2),  A(2),  ...,C(N), 
A(N),  D(l),  B(l),  D(2),  B(2),  ....  D(M),  B(M).  If  LO  =  1,  then  data  is  stored  sequentially  in  the 
following  variables:  N,  M,  LO,  AO,  BO,  P8,  A(l),  A(2),  ....  A(N),  B(l),  B(2),  ....  B(M);  the  arrays 
C(i)  and  D(j)  are  stored  with  ones  by  MLEQRE  in  this  case. 


AO  and  BO  contain  initial  estimates  of  A  and  B  (see  (82)),  supplied  by  the  user.  If,  however, 
BO  contains  zero  then  MLEQRE  supplies  the  initial  estimates. 


The  plotting  of  the  confidence  ellipses  at  the  50  and  95%  levels 
HP-9845  and  at  line  830  for  the  HP-85.  If  P8  >  2,  then  plot  appears  oi 
If  P8  =  1„  then  plot  appears  only  on  the  CRT.  If  P8  =  0,  then  no  plot  is  ci 


begins  at  line  1440  for  the 
CRT  and  also  the  printer, 
bnstructed. 


The  following  output  is  given  with  the  format  differing  slightly 
HP-85. 


between  the 


HP-9845  and 


Values  in  N,  M,  L0,  initial  values  in  A0,  B0,  P8. 

Values  of  c,,  a*  and  values  of  dj,  bj  if  LO  #  1 ,  otherwise  values  of  a{  and  bj  only. 

The  maximum  likelihood  estimates  MU(=#i),  SIG(=»o),  covariance  matrix  elements, 
initial  values  for  M=n/o),  B(=  1/a). 

Final  values  of  A0,  B0. 

Final  values  of  D1 ,  D2  ((N-R)  corrections). 

Number  of  (N-R)  inerations. 


X,  x  .  •*.  •*-.  .*  .*  +  x,  w.  **.  *•. 


/.  . 
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In  [7],  necessary  and  sufficient  conditions  were  derived  for  the  first  time  that  assure  the 
existence  and  moreover  the  uniqueness  of  the  maximum  likelihood  estimates  n,  a.  They  demand 
that 


max  bj  >  min  %  , 
j  i 


V.  1=1  j=i 


If  either  of  these  inequalities  is  not  satisfied,  a  message  is  printed  and  MLEQRE  terminates -for 
the  given  stimuli,  p  and  a  do  not  exist. 

A  Fortran  IV  program  based  on  [7]  and  upon  which  MLEQRE  is  modeled  has  been  available 
for  sometime  on  the  CDC  6700.  A  diluted  version  of  that  program  is  available  in  BASIC  for  the 
405 1-4054  series  Textronix  desk-top  computers,  [  13] .  It  does  not  contain  the  plotting  feature  for 
confidence  ellipses,  and  it  does  not  have  the  capability  to  take  advantage  of  the  increased  efficiency 
when  stimuli  are  repeated,  and  in  general  it  does  not  appear  to  be  as  efficient  as  MLEQRE.  The 
program  is  very  difficult  to  follow  and  we  have  not  been  able  to  verify  its  correctness. 


III.  LISTINGS  OF  PROGRAMS  AND  SAMPLE  OUTPUTS 

A  short  summary  of  the  input  and  output  associated  with  a  particular  program,  including 
examples,  is  given  starting  with  CIRCV.  This  is  followed  by  the  HP-9845  listing  for  that  particular 
program  and  then  the  corresponding  HP-85  listing. 

It  is  assumed  in  operating  the  HP-9845  programs  that 


PRINTER  IS  0 


has  been  executed. 


The  BASIC  language  for  the  HP-85  allows  for  multi-statement  lines,  where  the  statements  on  a 
numbered  line  are  separated  by  the  symbol  @.  Care  must  be  taken  in  interpreting  a  multi-statement 
line  when  an  IF  THEN  or  an  IF  •••  THEN  •••  ELSE  statement  is  not  the  last  statement  in  the 
line.  In  the  first  case,  when  the  IF-statement  is  true  and  THEN  is  not  followed  by  a  program 
transfer,  such  as  a  GO  TO,  the  execution  of  the  IF-statement  is  followed  by  executing  the  next 
statement  of  the  same  line.  When  the  IF-statement  is  false,  the  program  proceeds  directly  to  the 
next  sequentially  numbered  line.  For  example: 

100:  IF  A  =  B  THEN  C  *  B  @  GO  TO  500 
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I 


V* 

V* 


p 

S' 


If  A  =  B,  then  at  line  100  C  is  set  to  B  and  execution  continues  at  line  500.  If  A  ^  B  at  line  100, 
then  execution  continues  at  line  110. 

In  the  case  of  the  IF  •••  THEN  ELSE  statement,  if  the  IF  part  is  true,  and  THEN  is  not 
followed  by  a  program  transfer,  then  the  program  proceeds  to  the  next  numbered  line.  If,  however, 
the  IF  part  of  the  IF  •••  THEN  •••  ELSE  statement  is  false,  then  execution  of  the  ELSE  part  of 
the  statement  is  followed  by  executing  the  next  statement  of  the  same  line.  For  example: 

100:  IF  A  *  B  THEN  C  =  B  ELSE  D  =  B  @  GO  TO  500 
110:--- 

If  A  =  B,  then  C  is  set  to  B  and  program  proceeds  to  line  1 10;  if  A  #  B,  then  D  is  set  to  B  and  pro¬ 
gram  proceeds  to  line  500. 

CIRCV  or  GCEF 

Input:  R,  D,  V 

CIRCV:  R  =  R/ox  ,  D  =  v'x2  +  y */<xx  ,  V  =  1 . 

GCEF:  R  =  K  a  R /a.. ,  D  =  c  =  oy/ox  <  1 ,  V  »  0. 

If  R  and/or  D  <  0,  then  P  set  to  (-1 ).  Input  unacceptable. 

If  D  >  1  for  GCEF,  then  P  set  to  (-!).  Input  unacceptable. 

Output: 

CIRCV:  P  «  P(R,  D),  S2  =  || 

GCEF:  P  =  F(K,  c) ,  S2  =  (D/K)  3F/3K ,  D^0 

=  y/irj2  3F/3K ,  D  =  0 

Accuracy:  6  decimal-digits  for  P  and  S2 


Case 

EXAMPLES 

Input 

Output 

R 

y/x2  +  y2 

°x 

°y 

R 

D 

V 

P 

S2 

<D 

3 

4 

2 

2 

1.5 

2 

1 

.209232218046 

.214447041618 

© 

6 

4 

2 

2 

3 

2 

1 

.785637894167 

.101082811001 

© 

6 

5 

2 

2 

3 

2.5 

1 

.623010408440 

.130887896597 

© 

5 

6 

2 

2 

2.5 

3 

1 

.246101694960 

.130887896597 

© 

2 

0 

2 

4 

.5 

.5 

0 

.215288716030 

.738059987125* 

© 

8 

0 

2 

1 

4 

.5 

0 

.999926141520 

3.90547542875  (~5) 

© 

2 

0 

1 

0 

2 

0 

0 

.954499736146 

.135335283238 

© 

2 

0 

1 

1 

2 

1 

0 

.864664716770 

.135335283237 

*ox  and  CTy  inteichanged  jo  that  oy/ox  <  1 . 

Output  values  above  are  for  the  HP-85.  The  corresponding  values  for  the  HP-9845B  differ  in  the  last 
one  to  three  digits. 
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CIRCV —HP  9845 


LIST 

100  !  THIS  PROGRAM  IS  CALLED  “CIRCV".  IT  SUPPLIES  TWO 

FUNCTIONS:  P(R,D>,  THE  CIRCULAR  COVERAGE  FUNCTION 
105  !  OR  F<i<,  C> ,  THE  GENERALIZED  CIRCULAR  ERROR  FUNCTION. 

110  !  THE  INPUT  IS  R,D,V,  WHERE  IF  V=0  THEN  K=R  AND  C=D.  IF  V#0 

THE  OUTPUT  IS  P*P<R,D>;  IF  V  =  0  THEN  THE  OUTPUT  IS  P*F<K,C> 
115  !  INPUT  R  OR  D  <0  NOT  PERMITTED.  ALSO  FOR  V*0,  ABSCD-. 5> >. 5 

NOT  ALLOWED.  IN  SUCH  CASES  P  SET  TO  -1. 

120  !  LET  Pr  DENOTE  THE  PARTIAL  DERIVATIVE  OF  P  WITH 

RESPECT  TO  R.  THEN  Pr*R*S2. 

125  !  LET  Fk  DENOTE  THE  PARTIAL  DERIVATIVE  OF  F<K,C> 

WITH  RESPECT  TO  X.  THEN  Fk»<K/C>*S2,  Ct0. 

130  !  IF  C*0  THEN  Fk*  SQR<2/PI >*S2.  S2  IS  AVAILABLE 

INTERNALLY. 

135  !  SOURCES:  MATH  OF  COMP  APRIL  1961 , PP169, 173  AND  OCT. 

1961,  PP  375,382.  NWL  REPORT  #1768,  JAN.  1962. 

140  !  IEEE  TRANS.  INFO.  TH.  APRIL  1965,  P.  312. 

145  i  PROGRAM  IS  SET  FOR  SIX  DECIMAL  DIGIT  ACCURACY. 

150  P9<3)*21. 3853322378 

155  P3<2>*1. 72227577039 

160  P9<1>*. 316652890658 

165  Q9<2)*18. 9522572415 

170  Q9<  0*7. 8437457083 

175  R9<5)*7. 3738883116 

180  R9<4>*6. 8650184849 

185  R9<3>*3. 0317993362 

190  R9<2>*. 56316961891 

195  R9< 1 )*4. 3187787405E-5 

200  S9<4>*7. 3739608908 

205  S9<3>«15. 18490819 

210  S9<2>-12. 795529509 

215  S9<1)*5. 3542167949 

220  Bl*. 0000005 
225  S2-0 

230  P*0 

235  B2*.  70*7 106781 137 

240  IF  <R>-0>  AND  <D>*0)  THCN  255 

245  P*-l 

250  RETURN 

255  IF  R*0  THEN  250 

260  IF  V*0  THEN  605  !  COMPUTE  GCEF. 

265  Al-R-D 

270  A*ABS<AO 

275  IF  A<5. 386773  THEN  295 

280  IF  At<0  THEN  250 

285  P«i 

290  RETURN 

295  T-R*D 

300  T3*. 5*R*R 

305  B* . 5*D*D 

310  N-0 
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315  IF  T>7  THEN  455 

320  T 1 *B2*T  - 1 

325  T2=T3*B 

330  S0*EXP ( -T3-B  > 

335  S1“EXP<-B> 

340  IF  TS'.0005  THEN  355 
345  S 1 *S 1 *T3 

350  GOTO  36l  , 

355  S 1 =S 1 -30 

360  S2“S0 

365  T0*S  1 

370  N*N+ 1 

375  M»ON 

380  S0“T2*M*M*S0 

385  T0*B*M*T0-S0 

390  S1=*S1+T0 

395  S2*S2+S8 

400  IF  TON  THEN  370 

405  IF  T0>B1  THEN  420 

410  P*S 1 

415  RETURN 

420  N-N+l 

425  W*l/N 

430  S0*T2*H*M*S 

435  T0«B*M*T0-S0 

440  S1-S1+T0 

445  S2-S2+S0 

450  GOTO  405 

455  Tl-2*flBS<T3-B> 

460  fl»fl*B2 
465  T3-1/CT+T) 

470  T2-S0RCT3) 

475  Sl».5**l*fll 
480  S2-EXP<-S1> 

485  S0-. 564189503545#T2#S2 

490  GOSUB  710 

495  T0»<R+D>*B2#T2*E 

500  T2*S1*T3 

505  T3-.5*T3 

510  S1-T0 

515  S2-S0 

520  N-N+2 

525  M»N- 1 

530  fl»M'N 

535  S0-fl«T3*S0 

540  T0»Tl*S0-T2*fl*T0 

545  S0»M*S0 

550  S1-S1+T0 

555  S2-S2+S0 

560  IF  T0-B1 >0  THEN  520 

565  IF  S0-B1 >0  THEN  580 

570  P«. 5#RBS< l+SGN<fll >-S2- :GN<fll >#S1 > 

575  RETURN 
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580  N=N+2 

585  M*N- 1 

590  S0*M*M*T3*S0/'H 

595  S2*S2+S0 

600  GOTO  565 

605  IF  flBS<D-. 5) >. 5  THEN  245  ! START  FOR  GCEF . 

610  IF  R>-5. 386773  THEN  285 

615  IF  D<  >0  THEN  650 

620  ft»B2*R 

625  Sl»fi*fi 

630  S2»EXP<-S1> 

635  GOSUB  710 
640  P*l-E 
645  RETURN 
650  K*R 
655  C*D 
660  T-.5xC 
665  R*K*<1+C>*T 

670  D*K*< 1-C)*T 

685  GOSUB  265 
690  R-K 
695  D-C 

700  P*flSS<P+P+S2-i) 

705  RETURN 

710  IF  flBS<fl>>.5  THEN  725  !E»ERFC<fl> 

715  E»l~fl*<<P9<l>*Sl*P9<2>  >#S1*P9<3>  >^<  CS1*Q9<1>  )*S1+Q9<2)  > 

720  RETURN 

725  E»<<<<R9<l>*fl+R9<2>>*fi+R9<3>>«fl*R9<4>>*ft+R9<5>>'<<<<fl+S9<l>>*«+S9<2> 

<3> )*ft+S9<4>  >*S2 
730  RETURN 
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100  !  THIS  PROGRAM  IS  CALLED  “Cl 
RCV" .  IT  SUPPLIES  TWO  FUHCTI 
ONS :  P<R,D>,  THE  CIRCULAR  CO 
VERAGE 

105  !  FUNCTION  OR  F<K,C>  THE  GEN 
ERALIZED  CIRCULAR  ERROR  FUNC 
TION. 

110  !  THE  INPUT  IS  R,D.V,  WHERE 
IF  V=0  THEN  K=R  AND  C=D.  IF 
V#0  THE  OUTPUT  IS  P=P<R,D>; 
IF  V=0 

115  !  THEN  THE  OUTPUT  IS  F<K,C>. 
INPUT  R  OR  D<8  NOT  PERMITTED 
ALSO  FOR  V=e  ABS<D- . 5)> .5  N 
OT 

120  !  ALLOWED.  IN  SUCH  CASES  P  S 
ET  TO  -1 . 

125  !  LET  Pr=THE  PARTIAL  DERI VAT 
IVE  OF  P  WITH  RESPECT  TO  R. 
THEN  Pr=R*S2. 

139  !  LET  Fk=  THE  PARTIAL  DERIVA 
TIVE  OF  F  WITH  RESPECT  TO  K. 
IF  C#8  THEN  Fk=<K/C)*S2 . 

135  !  IF  C=0  THEN  Fk=SQR<2^PI >*S 
2.  S2  lb  HVHILABLE  INTERNALL  Y 
148  !  PROGRAM  IS  SET  FOR  6-DIGIT 
ACCURACY . 

145  '  SOURCES :  MATH  OF  COMP.  APR 
IL, 1961 iPP. 169—173  AND  OCT. 
1961,  PP.  375-382. 

150  !  SOURCES ;  NWL  REP0Rt#1768, J 
AN. 1962.  IEEE  TRANS.  INFO.  T 

H.  APRIL, 1965,  P.312. 

155  P9<3)=21 .3853322378  8  P9<2>« 

I. 72227577039  8  P9<1>=  31665 
28^0658 

160  Q9<2>=18. 9522572415  8  Q9<1>= 
7  8437457083 

165  R9<5>=7. 37388831 16  8  R9<4)=6 
.8650184849  8  R9<3>=3 . 021799 
3362  8  R9<2)=. 56316961891 
178  R.9<  1  >=4 . 3 1 87787405E-5  8  S9<4 
) =7. 3739608908  8  S9<3)*»5.18 
490819  8  S9<2)=12  795529509 
175  S9<1>=5  3542167949 
188  Bl=  8808005  8  S2=0  8  P=0  8  E 
=88  B2=  707186781187 
185  IF  R>=8  AND  D>=0  THEN  195 
190  p-_!  0  RETURN 
195  IF  R=0  THEN  RETURN 
208  IF  V=0  THEN  325 
205  A1=R-:  ••  A=ABS < A 1 )  8  IF  A<5. 

386773  T HEN  220 
210  IF  A 1 <0  THEN  P=0  ELSE  P=1 
215  RETURN 

220  T=R*D  8  ^7=  5*r*r  @  B=  5*0*0 
8  N=0 


2.25  IF  T> 7  THEN  275  ELSE  T1^B2*T 
-1 

230  T2=T3*B  8  S0=EXP<-T3-B>  8  IF 
T3> . 0095  THEN  SI =EXP  .  -B)-SO 
ELSE  S1=EXP<-B>*T3 
235  S2=S9  8  T0=S1 
248  N=N+1  8  M=l/N 
245  S0=T2*M*M*S0  8  T0=B*M*T0-S0 
258  S1=T0+S1  8  S2=S2+S8 
255  IF  T1>N  THEN  248 
260  IF  T0>B1  THEN  278 
263  P=cl  8  RETURN 
270  N=N+1  8  M=l/N  8  S0=T2*M*M*S0 
8  T8=B*H*T0-S0  8  S1=T0+S1  8 
S2=S2+S6  8  GOTO  260 
275  T1=2*ABS<T3-B>  8  A=A*B2  8  T7 
=1'<T+T>  8  T2=SQR<T3>  8  Sl=. 
5*A1*A1 

280  S2=EXP(-S1>  8  S0=.5*T2*1  128 
37916709*S2  8  T0=<R+D>*B2*FN 
E<A)*T2 

2b5  T2=S1*T3  8  T3=.5*T3  e  S1=T0 
8  S2=S8 

299  N=H+2  8  M=N-1  8  A=M'N 
295  S8-A*T3*S0  8  T0=T1*S8-T2*A*T 
0  8  S0=M*S0 

308  S1=S1+T0  8  S2=S2+S8 
305  IF  T0— B1>0  THEN  298 
318  IF  S0-B 1 > 0  THEN  320 
315  P= . 5* ABS < 1 +SGN  <  A 1 > -S2-SGN <  A 1 
>*S1)  8  RETURN 

320  N=N+2  8  M=N-1  8  S8=M*M*T3*S0 
y'H  8  S2-S2+S0  8  GOTO  310 
325  IF  ABS<D-.5><=.5  THEN  335 
339  P=-l  e  RETURN 
335  IF  R <5  386773  THEN  345 
346  P=1  8  RETURN 
345  IF  Dt8  THEM  355 
350  A=B2*R  e  S1=A*A  8  S2=EXP<-S1 

>  8  P*1-FNE<A>  8  RETURN 

355  K=R  8  C=D  8  T*.5/'C  8  R=K*<1+ 
C)*T  8  0»K*<1-C>*T 
368  GOSUB  205 

365  R=K  e  0=C  8  P*ABS<P+P+S2-i > 
370  RETURN 

375  DEF  FNE< A>  !  FNE < A > *ERFC < A > 
380  IF  ABS<A>> .5  THEN  395 
385  FNE=1-A*<<P9< 1 )*S1+P9<2> )*S1 
♦P9<3>  >/<  <S1+Q9<  1  >  >*Si+Q9<2) 

> 

390  GOTO  480 

395  FNE*< <  <  <R9< 1 >*A+R9<2>  >*A+R9< 
3>  >*A+R9<4)  )*A+R9<5>  >/'<<<  (A+ 
S9< 1 >)*A+S9<2) >*A+S9<3) )*A+S 
9<4) )*S2 
480  FN  END 
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Input:  R,  H,  R,  SI ,  S2 

(H  =  h,  K  =  k,  SI  =  ox,  S2  =  oy) 

Output:  P  =  P(R,  h,  k,  ox,  oy) 

Accuracy:  6-decimal-digits  for  ELLCV 
3-decimal-digits  for  ELLCV3 


EXAMPLES 

1)  R  =  5,  H  *  2,  R  =  3,  ox  =  SI  =  3,  oy  =  S2  =  2  (See  Example  3  of  POLYCV,  page  59) 

P  =  .588490575749  (ELLCV) 

P  =  . 5 8849057921 7  (ELLCV3)  ■ 

2)  R  =  2,  fl  =  0,  K  =  0,  ax  =  2,  oy  =  4  (See  Example  ©  for  CJRCV,  page  26) 

P  =  .2 1 5 2887 1 6038  (ELLCV) 

P  =  .215288734652  (ELLCV3) 

3)  R  =  3,  H  =  2.  R=  2  y/7,  ax  =  2,  oy  =  2  (See  Example  ©  for  CIRCV,  page  26) 

P  =  . 209232220601  (ELLCV) 

P  =  .209232478927  (ELLCV3) 

4)  R  =  3,  fl  =  0,  R  =  .3,  ox  =  1,  oy  =  1/10 

P  -  .997146500681  (ELLCV) 

P  =  .99709445 1746  (ELLCV3) 


NOTE:  The  bare  on  R,  H,  K  are  to  conform  with  the  discuaion  in  the  previous  section,  As  BASIC  variables  the 
bars  are  deleted. 
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100  !  THIS  PROGRAM  IS  CALLED  "ELLCV".  IT  SUPPLIES 

THE  ELLIPTICAL  COVERAGE  FUHCT I  ON : P < R , H, K,  S 1 ,  S2 >  . 

105  !  P  DENOTES  THE  PROBABILITY  OF  A  SHOT,  NORMALLY 

DISTRIBUTED  WITH  MEAN  <0,0>  AND  STANDARD 
110  !  DEVIATIONS  S1,S2  IN  THE  X  AND  Y  DIRECTIONS, 

RESPECTIVELY,  FALLING  IN  A  CIRCLE  IN  THE 
115  !  XY-PLANE  OF  RRDIUS  R  AND  CENTERED  AT  <H,K>.  THE 

INPUT  IS  R,H,K,S1,S2.  THE  OUTPUT  IS  P. 

120  !  PROGRAM  IS  SET  FOR  6-DECIMAL-DIGIT  ACCURACY  IN  P. 

125  !  ELLCV  USES  ERFC  WITH  9  DIGIT  RELATIVE  ACCURACY. 

130  !  SOURCES:  NWL  REPORT  #1710,  AUG. I960.  MATH  OF 

COMP.  OCT.  1961,  PP.  375,382. 

135  !  INPUT  R,H,K,S1,S2 

140  !  "ELLCV"  CONSTRUCTED  IN  COLLABORATION  WITH  ALFRED  MORRIS. 

145  OPTION  BASE  1 

150  DIM  P9<3),Q9<2)jR9<5),S9<4),X(43),Y<43> 

155  P9< 1 )*3. 1 6652890658E- 1 

160  P9<2)«1. 72227577039 

165  P9<3)*21. 3853322378 

170  Q9< 1)*7. 8437457083 

175  Q9<2>»18. 9522572415 

180  R9< 1 )*4. 3197787405E-5 

185  R9<2!>». 56316961891 

190  R9<3>-3. 0317993362 

195  R9<4>*6. 8650184849 

200  R9<5>-7. 3738883116 

205  S9(l)»5. 3542167949 

210  S9<2>-12. 795529509 

215  S9<3>«15. 18490819 

220  S9< 4 >■ 7. 3739608908 

225  XU)-.  238619186083 

230  X<2>». 661209386466 

235  X<3>-. 932469514203 

240  X<4>». 183434642496 

245  X<5>*. 525532409916 

250  X<6>«. 796666477414 

255  X<7>«. 960289856498 

260  X<8)». 125233408511 

265  X<9)-. 367831498998 

270  X<10)-. 587317954287 

275  Xai>«.  769902674194 
280  X<12>-. 90411725637 

285  X< 13>-. 981560634247 
290  X<14>-9.50125898376E-2 

295  X< 15)». 281603550779 
300  X(16)». 458016777657 

305  X<17)-. 617876244403 

310  X< 18>». 755404408355 
315  X<19)-. 865631232388 

320  X<20>». 944575023073 

325  X<21>». 909400934992 
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330  X<22)*7.65265211335E-2 

335  X<23)=. 227785851142 

340  X<24)*. 373706088715 

345  X<25)*. 510867801951 

350  X<26>“3 636053680727 

355  X<27)*. 74633190646 

360  XC28)». 839116971822 

365  X<29)*=. 912234428251 

370  X< 30)*. 963971927278 
375  X<31)-. 993128599185 

380  X(32)*6. 4056892 86 26E -2 

385  X<33)-. 191 1 18867474 

390  X<34)*. 315042679696 

395  X<35>*. 433793507626 

400  X<36)-. 545421471389 

405  X<37)«. 648893651937 

410  X<38)». 740124191579 

415  X<39)», 820001985974 

420  X<40)-. 886415527004 

425  X(41>*. 938274552003 

430  X<42>*. 974728555971 

435  X<43>-. 995187219997 

440  Y<1>«. 467913934573 

445  Y<2)». 360761573048 

450  Y<3>-., 171324492379 

455  Y<4>*. 362683783378 

460  Y<5)«. 313706645878 

465  Y<6)-. 222381034453 

470  Y<7)-. 10122853629 

475  Y<8)-. 249147045813 

480  Y<9)-. 233492536538 

485  Y< 10)-. 203167426723 
490  Y(ll)., 160078328543 

495  Y<12>*. 106939325995 

500  Y< 13>-4. 7 1 753363865E-2 
505  Y<14>*. 189450610455 

510  Y< 15)-. 182603415045 
515  Y<16>*. 169156519395 

520  Y<1?>*. 149595988817 

525  Y< 18)-. 124628971256 
530  Y<19>«9. 515851 16825E -2 

535  Y<20)-6. 22535239386E-2 

540  Y<21)-2.71524594U8E-2 

545  Y<22>-. 152753387131 

550  Y<23>-. 149172986473 

555  Y<2*>". 142096109318 

560  Y<25)-. 131688638449 

565  Y(26>«. 118194531962 

570  Y<27>-. 101930119817 

575  Y< 28) -8. 3276741576 7E -2 

500  Y<29>-6. 2672048334. E-2 
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£35  Y  ■>'  39  >  =  4  . 36014298004E-2 

593  Y<31>=1.761400?1392E-2 

595  v(32)=. 127938195347 

600  Y< 33 >=.12583745634? 

605  7<34>». 121678472928 

610  Y<35>=. 115505668054 

615  Y<36>=. 107444270116 

620  Y<3?'»=9.?6186521041E-2 

625  Y(38)=. 086190161532 

630  Yt39>*7.33464814111E-2 

635  Y<40>-5. 92985 049154E-2 

640  Y<41)=4.42??4388174E-2 

645  Y<42)-2.85313886289E-2 

650  Y<43>«. 0123412298 

655  8=4.892 

660  81*5.387 

665  82=3.8775 

670  83=5.16 

575  Bl=. 564189583548  !  1/S3R<PI> 

680  B=l. 41421356237  !  SQR<2) 

685  ,  B2-29.0197C9  !B2-81*81 

690  P=0 

695  23=.000001*S1*S2 

700  IF  R*R<  =23  THEN  RETURN 

705  H2-H*H+K*K 

710  D-M8X ( S 1 , S2  > 

715  T*R-81*D 

720  !  PROCEED  TO  SEE  IF  P-0  OR  P-1 

725  IF  T< 0  THEN  745 

730  IF  T*T<H2  THEN  745 

735  P-1 

740  RETURN 

745  H8-8BS  <  H ) 

750  K8-8BSOO 

755  IF  R-H8+8*S 1 <  =0  THEN  RETURN 

760  IF  R-K8+8*S2< -0  THEN  RETURN 

765  S0-S3R<H2> 

770  IF  S10S2  THEN  790 

775  H8-S0 

780  K8-0 

785  IF  R+8*S  1  '<  -HO  THEN  RETURN 

790  IF  S0< -R  THEN  910 

795  D-CSO-R^D 

800  IF  R#R*EXP<-.5*D#D>>23  THEN  910 
805  RETURN 

810  IF  S0<R+81*MIN<S1 , S2>  THEN  910 

815  IF  H8*K8-0  THEN  910 

820  H9-H8/S1 

825  K9-K8/S2 

830  D-H9*H9+K9*K9 

835  IF  D< -B2  THEN  910 
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840  Z2-R/S2 

845  Q-S2'S1 

850  Q1*C*Q 

855  F*Q1*H9*H9+K9*K9 

860  Z1-22*Z2*F/-D 

865  Z-D-Z1-B2 

870  IF  2<0  THEN  880 

875  IF  2*2-4*21 *B2>*0  THEN  RETURN 

880  T1*H8*H8+Q1 *K8*K8 

885  28»B2*S1*S1*T1/H2 

890  R2«R*R 

895  Y-H2-R2-Z8 

900  IF  Y<-0  THEN  910 

905  IF  Y*Y-4*R2*Z8>«0  THEN  RETURN 

910  28*0!  FINE  LIMITS  OF  INTEGRATION 

915  Z*K8+A3*S2 

920  H3*K8-A3*S2 

925  S0*S1 

930  S9*S2 

935  2*R-2 

940  D 1 *0 

945  IF  2>*0  THEN  D1*SQR<2/R> 

950  IF  H3>*0  THEN  970 

955  H5*0 

968  E3-1-D1 

965  GOTO  980 

970  E3*SQR<1-H3/R)-D1 

975  H5-1 

980  IF  Z8O0  THEN  1020 

985  28-1 

990  F-E3 

995  T-D 1 

1000  Z-H8+A3*S1 

1005  H6-H5 

1010  H3*H8-A3*S 1 

1015  GOTO  935 

1020  IF  F>-F.3  THEN  1065 

1025  E3-F 

1030  Dl-T 

1035  S9-S1 

1040  Z8-H8 

1045  S0-S2 

1050  H8-K8 

1055  K8-28 

1060  H5-H6 

1065  E3*.5*E3  !  BEGIN  GAUSSIAN  INTEGRATION 

1070  N-E3*«*<.34/S0+l/'<.025*ABS(R-’<8>+5*S9)) 

1075  Z2-R/<B*S9> 

ie80  R8*R/<B*S0> 

1085  H0*H8/,<B*S0) 

1090  K8*K8»'<B*S9> 
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1095  IF  N<2. 75  THEN  1115 

1100  J-31 

1105  N 1  ■  1 2 

1110  GOTO  1205 

1115  IF  NC1.35  THEN  1135 

1120  J-21 

1125  Ni-10 

1130  GOTO  1205 

1135  IF  N< . 75  THEN  1155 

1140  J-13 

1145  Nl-8 

1150  GOTO  1205 

1155  IF  N< . 35  THEN  1175 

1160  J-7 

1165  Nl-6 

1170  GOTO  1205 

1175  IF  N< . 15  THEN  1195 

1180  J-3 

1185  Nl-4 

1190  GOTO  1205 

1195  J-0 

1200  Nl-3 

1205  2-Z3-0 

1210  Y«B1*E3*R8 

1215  K9-1 . 04E-8 

1220  H9-1. 9999999792 

1225  G3-0 

1230  M-N1+N1 

1235  I— N1 

1240  IF  K8-0  THEN  1415 

1245  FOR  L"  1  TO  H 

1250  IF  1-0  THEN  1-1 

1255  T-E3*:SGN< I >#X< J+RBS< I >  >  +  l >+Dl 

1260  T9-T  *T 

1265  T 1 -R3*< 1 -T9 ) 

1270  T2-H8-T 1 

1275  T4-EXP<-T2*T2> 

1280  IF  H8O0  THEN  1295 

1285  T4-T4+T4 

1290  GOTO  1310 

1295  IF  H5O0  THEN  1310 

1300  T2«H8fTl 

1305  T4-T4+EXP<rT2#T2> 

1310  IF  ZO0  THEN  1340 

1315  Z1-Z2#T*SQR<2-T9> 

1320  K1-K8-Z1 

1325  IF  RBS(K1  XR2  THEN  1350 

1330  IF  K 1 >0  THEN  1395 

1335  Z-l 

1340  K5-H9 

1345  GOTO  1360 
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1350  GOSUB  1560 

1355  K5=K3 

1360  K 1 =K8+Z 1 

1365  IF  K 1  <  H2  THEN  1380 

1370  K5=K5-K9 

1375  GOTO  1390 

1380  GOSUB  1560 

1385  K5*K5-K3 

1390  G3=G3+K5*T4*T*Y<J+HBS<I>> 

1395  1=1+1 

1400  NEXT  L 

1405  P=Y*G3 

1410  RETURN 

1415  FOR  L* 1  TO  H 

1420  ,  IF  1*0  THEN  1=1 

1425  T=E3*<SGN<I>*X<J+flBS<I>>+l>+Dl 

1430  T9=T*T 

1435  T1»R8*< 1-T9> 

1440  T2-H8-T1 

1445  T4=EXP<-T2*T2> 

1450  IF  H8O0  THEN  1465 

1455  T4-T4+T4 

1460  GOTO  1480 

1465  ,  IF  H5O0  THEN  1480 

1470  T2*H8+T 1 

1475  T4«T4+EXP<-T2*T2) 

1480  IF  Z<  >0  THEN  1500 

1485  K1»Z2*T*SQR<2-T9> 

1490  IF  Kl<02  THEN  1510 

1495  Z=1 

1500  K5=H9 

1505  GOTO  1520 

1510  GOSUB  1560 

1515  K5=2*  < 1 -K3  > 

1520  G3=G3+K5*T4*T»Y< J+flBS( I >  > 

1525  1=1+1 

1533  NEXT  L 

1535  P=Y»G3 

1540  RETURN 

1545  REN  CODY  FOR  K3=ERFC<K1 >-~9  DIGITS 

1550  !  IF  Kl<—  02  THEN  K3-2-2.08E-8 

1555  !  IF  K 1 >=82  THEN  K3=l,04E-8 

1560  IF  RBS<K1 )  >.  5  THEN  1580 

1565  K4-K1+K1 

1570  K3=1-K1*<<P9<1  )+K4+P9<2>  )*K4+P9<3>  )>■'<  <K4+Q9<  1  >  >*K4+Q9<2)  > 

1575  RETURN 

1580  K4=flBS<Kl) 

1590  X3»><C<<R9<1)*K4+R9<2) >#K4+R9<3> >*K4+R9<4> >*K4+R9<5> )/<  <  < (K4+S9 

<1>>»K4+S9<2>>*K4+S9<3>>#K4+S9<4>>#EXP<-K1*K1> 
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168a  RETURN 

1605  !  K6«1/'<K1#K1).  NOT  USED  PRESENTLY— FOR  USE  WHEN  Kl>4. 

1610  !  K3»<B1+K6*<<<V9<1)#K6+V9<2))*K6+V9<3))*K6+V9<4>)/<<<K6+W9<1))* 

K6+W9<2>>*K6*W9<3>>>*EXP<-K1*K1 >/K4 

1620  X7-0  !  CRUTCHER  TABLE  CHECK 

1625  INPUT  R,H,K,S1,S2 

1630  COSUB  100 

1635  print  r;k;k;si;S2 

1640  R4»R 

1645  X7-X7+.1 

1650  R-R4*X7 

1655  GOSUB  690 

1660  INACE  DD. DD, 2X, D. DDDDDD 

1665  PRINT  USING  1660JX7;P 

1670  IF  P<«. 999999  THEN  1645 

1673  BEEP 

1680  END 
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100  !  THIS  PROGRAM  IS  CALLED  “  ELLCV3" .  IT  SUPPLIES 

THE  ELLIPTICAL  COVERAGE  FUNCT I  ON : P < R , H , K , S 1 , S2 > . 

105  !  P  DENOTES  THE  PROBABILITY  OF  A  SHOT,  NORMALLY 

DISTRIBUTED  WITH  MEAN  <0,0>  AND  STANDARD 
110  !  DEVIATIONS  S1,S2  IN  THE  X  AND  Y  DIRECTIONS, 

RESPECTIVELY,  FALLING  IN  A  CIRCLE  IM  THE 
115  !  XY-PLANE  OF  RADIUS  R  AND  CENTERED  AT  CH,K>.  THE 

INPUT  IS  R,H,K,Si,S2.  THE  OUTPUT  IS  P. 

120  !  PROGRAM  IS  SET  FOR  3-DECIMAL-DIGIT  ACCURACY  IN  P. 

125  !  ELLCV3  USES  ERFC  WITH  9  DIGIT  RELATIVE  ACCURACY. 

130  !  SOURCES:  NWL  REPORT  #1710,  AUG. I960.  MATH  OF 

COMP.  OCT..  1961,  PP.  375,382. 

135  !  INPUT  R,H,K,S1,S2 

140  !  “ELLCV3"  CONSTRUCTED  IN  COLLABORATION  WITH  ALFRED  MORRIS. 

145  OPTION  BASE  1 

150  DIM  P9<3),Q9<2),R9<5),S9<4),X<21>,Y(21> 

155  P9<1>*. 316652890658 

160  P9<2>»1. 72227577039 

165  P9<3>*21. 3853322378 

170  Q9<1W.  8437457083 

175  Q9<2>-18. 9522572415 

180  R9<l)*4.3187787405E-5 

185  R9 <2 >*.56316961891 

190  R9<3>«3. 0317993362 

195  R9<4>«6. 8650184849 

200  R9<5>*7. 3738883116 

205  S9<1>*5. 3542167949 

210  S9<2>*12. 795529509 

215  S9<3>*15. 18490319 

220  S9<4>-7. 3739608908 

225  X<1 >*. 238619186083 

230  X<2>*. 661209386466 

235  X<3>*. 932469514203 

240  X<4>*. 183434642496 

245  X<5>*. 52553240991,6 

250  X<6>*. 796666477414 

255  X<?>*. 960289856498 

260  X<8>*. 125233408511 

265  X<9>*. 367831498998 

270  X< 10>*. 587317954287 
275  X  < 1 1 > - . 769902674194 

280  X<12>“. 9041172563? 

285  X<13>«. 981560634247 

290  X<14>»9.50125098376E-2 

295  X<15>*. 281603550779 

300  X<16>*. 458016777657 

305  X<17>«. 617876244403 

310  X<18>*. 755404408355 

315  X<19>*. 865631202388 

320  X<20>*. 944575023073 

325  X<21>*. 909400934992 

330  Y<1>*. 467913934573 

335  Y<2>*. 360761573048 

340  Y<3>*. 171324492379 
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345  Y<41-. 362683783378 

350  Y<51=. 313706645878 

355  YC61-. 222381034453 

360  YC71-. 10122853629 

365  Y<8>». 249147045813 

370  Y<9>-. 233492536538 

375  Y<10>-. 203167426723 

380  Y<11>». 160878328543 

385  Y<12>-. 106939325995 

390  Y<131-4.71753363865E-2 

395  Y  <  1 4  >  - . 189450610455 

400  Y< 131-. 182603415045 
405  Y<16>». 169156519395 

410  YC 171*. 149595988817 
415  Y<18>*. 124628971256 

420  Y<19>*9. 515851 16825E-2 

425  YC201-6.22535239386E-2 

430  Y<21)*2. 71524594 118E-2 


435 

fl-3.291 

440 

Al-3. 89895 

443 

A2-2.898 

450 

A3-3.7 

453 

B1-. 564189583548 

!  1/SQR<PI  > 

460 

B-l. 41421356237 

!  SQR  <21 

465 

B2-13. 2018111 

! B2-R 1 1 

470 

P-0 

473 

Z3«.001*S1#S2 

480  IF  R*R<»Z3  THEN  RETURN 
485  H2»H#H+K*K 
490  D«MAX<S1,S2> 

495  T-R-AHU 

500  !  PROCEED  TO  SEF.  IF  P-0  OR  P-1 

505  IF  T<0  THEN  525 

510  IF  T*T<H2  THEN  525 

515  P-1 

520  RETURN 

525  H8-ABSCH1 

530  K8»ABS<IO 

535  IF  R-H8+A*Sl<«0  THEN  RETURN 

540  IF  R-K8+A*S2<«0  THEN  RETURN 

545  S0-SQR<H2) 

550  IF  S10S2  THEN  570 
553  H8-S0 

560  K8-0 

565  IF  R-H8+A*Sl<»0  THEN  RETURN 

570  IF  S0<-R  THEN  690 

373  D-<S0-R>/-D 

580  IF  R*R*EXP<-.3#D#D>>Z3  THEN  690 
583  RETURN 

590  IF  S0<R+A1#MIN<S1,S2>  THEN  690 

595  IF  H8#K8-0  THEN  690 

600  H9-H8/SI 
603  (O-KS/'SZ 
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610  D-H9*H9+K9*K9 

615  IF  D<»B2  THEN  690 

620  Z2-R/S2 

625  Q-S2/S1 

630  Q1-Q*Q 

635  F-Q1*H9*H9+K9*K9 

640  Z1-Z2*Z2*F/D 

645  Z-D-Z1-B2 

650  IF  Z<0  THEN  660 

655  IF  Z#Z-4*Z1*B2>»0  THEN  RETURN 

660  T1«H8*M84-Q1*K8*K8 

665  Z8«B2*S1*S1*T1^H2 

670  R2-R#R 

675  V-H2-R2-Z8 

680  IF  Y<-0  THEN  690 

685  IF  Y*Y-4*R2*Z8>*0  THEN  RETURN 

690  Z8-0  !  FIND  LIMITS  OF  INTEGRATION 

695  Z*K8+R3»S2 

700  H3»K8-N3*S2 

70S  S0-S1 

710  S9-S2 

715  Z-R-Z 

720  Dl-0 

725  IF  Z>0  THEN  D1'SQR<Z/R> 

730  IF  H3>-0  THEN  750 

735  H5-0 

740  E3-1-D1 

745  COTO  760 

750  E3-SQR<1-H3<'R)-D1 

755  H5-1 

760  IF  Z8O0  THEN  800 

765  Z8-1 

770  F-E3 

775  T«Dl 

780  Z*H8*R3*S1 

785  H6-H5 

790  H3-H8-A3#S1 

795  GOTH  715 

800  IF  F>»E3  THEN  845 

805  E3-F 

810  Dl-T 

815  S9-S1 

820  Z8-H8 

825  S0-S2 

830  H8-K8 

835  K8-Z8 

840  N5-H6 

845  E3-.S*E3  !  BEGIN  GAUSSIAN  INTEGRATION 

850  N«E3*R*<.34^S0<M/<  . 025#ABS<R-K8>*5*S9>  > 

855  Z2-R^<B*S9> 

860  R8-R/<B*S8> 

865  H8«N8/’<B»S0> 

870  K8-K8s<B»S9) 
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873  IF  N<2  THEN  895 

880  J-13 

883  Nl-8 

890  GOTO  943 

895  IF  H< . 675  THEN  913 

900  J-7 

905  Nl-6 

910  GOTO  945 

915  IF  N< . 3  THEN  933 

°20  J-3 

923  Nl-4 

930  GOTO  945 

933  J-0 

940  Nl-3 

945  Z-Z3-0 

930  Y«B1*E3*R8 

953  K9-. 000013 

960  H9«l. 99997 

965  G3-0 

970  N-Nl+Nl 

975  I— HI 

980  IF  K8-0  THEN  1160 

985  Z3-0 

990  FOR  L»1  TO  H 

993  IF  1-0  THEN  I-J 

1000  T«E3*  <  SGN< I )*X< I >X1  >+Dl 

1005  T9-T*T 

1010  Tl-R8#< 1-T9) 

1013  T2-HQ-T 1 

1020  T4»EXP<-T2#T2) 

1023  IF  H8<  >0  THEN  1040 

1030  T4-T4+T4 

1033  GOTO  1033 

1040  IF  H3<  >0  THEN  1033 

1043  T2-H8+T 1 

1030  T4-T4*EXP<-T2»T2> 

1033  IF  ZO0  THEN  1083 

i860  Z1-Z2*T#SQR<2-T9> 

1063  Kl-KS-Zl 

1070  IF  PBS<KlXfl2  THEN  1093 

1073  IF  K1 >0  THEN  1100 

1080  Z-l 

1083  K3-H9 

1090  GOTO  1105 

1093  GOSUB  1300 

1100  K3-K3 

1103  K1-K8+Z1 

1110  IF  Kl<fl2  THEN  1123 

1113  K3-K3-K9 

1120  GOTO  1133 

1123  GOSUB  1303 

1130  K3-K3-K3 

1133  G3-G3^K3*T4#T#Y< J*flBS< I ) ) 
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1148  I-I+l 

1143  NEXT  L 

1130  P«Y#G3 

1135  RETURN 

1150  FOR  L«1  TO  H 

1163  IF  1-0  THEN  1-1 

1170  T»F3*<SGN<I>*X<J*RBS<I>>+1>*D1 

1173  T9«T*T 

1180  T1»R8*<1-T9> 

1183  T2-H8-T1 

1190  T4»£XP<-T2*T2> 

1193  IF  H8<  >0  THEN  1210 

1200  T4-T4+T4 

1203  GOTO  1225 

1210  IF  H3O0  THEN  1223 

1215  T2-H8+T1 

1220  T4-T4+EXP<-T2«T2> 

1225  IF  ZOO  THEN  1245 

1230  K1»Z2*T*SQR<2-T9> 

1235  IF  K1<R2  THEN  1233 

1240  Z-l 

1243  K3-H9 

1230  GOTO  1263 

1233  GOSUB  1303 

1260  K3-2*<l-K3> 

1263  G3**G3+K3*T4*T*Y<  J*RBS<  I  >) 

1270  I*I+1 

1273  NEXT  L 

1280  P«Y#G3 

1283  RETURN 

1290  REN  CODY  FOR  K3-ERFC<K1  >—  9  DIGITS 
1293  ,!  IF  Kl<— R2  THEN  K3-H9. 

1300  I  IF  K1 >»R2  THEN  K3»K9. 

1303  IF  R8S<K1>>.5  THEN  1323 

1310  K4-K1#K1 

1313  K3»1-K1*<  <P9<  1>*K4*P9<2>  >#K4«-P9<3>  J'C  <K4+Q9<  1  >  >*K4+Q9<2)  > 

1320  RETURN 

1323  K4«RBS<K1> 

1330  K3«<<  <  <R9< U#K4*R9<2>  >*k4+R9<3> >*K4+R9<4>  >#K4«-R9<3> >/<  < <  <K*+S9 

< l > >#K4*S9<2> >*K4*S9<3>  >*K4*S9<  4>  >*EXP< -Ki »K1 > 

1333  IF  K 1 < 0  THEN  K3-2-K3 

1340  RETURN 

1343  X7«0  !  CRUTCHER  TRBLE  CHECK 

1330  INPUT  R,H,K,S1,S2 

1333  GOSUB  100 

1360  PRINT  r;h;k;si;S2 

1363  R4-R 

1370  X7»X7*. 1 

1373  R-R4*X7 

1380  GOSUB  478 

1383  INRGE  DD. DD, 2X, D. DDDD 

1390  PRINT  USING  1383;X7;P 

1393  IF  P<-. 99999  THEN  1370 

1400  BEEP 
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lee  !  THIS  PROGRAM  IS  CALLED"ELL 
CVMT  SUPPLIES  THE  ELLIPTIC 
AL  COVERAGE  FUNCTION ’ P<R,H,K 
,S1,S2> 

195  !  P  DENOTES  THE  PROBABILITY 
OF  A  SHOT ,  NORMALLY  DISTRIBU 
TED  WITH  MEAN  <0,0>  AND  STAN 
OARD 

110  »  DEVIATIONS  SI, 32  IN  THE  X 
AND  Y  DIRECTIONS,  RESPECTIVE 
LY,  FALLING  IN  A  CIRCLE  IN  T 
HE 

115  !  XY-PLANE  OF  RADIUS  R  AND  C 
ENTERED  AT  <H,K>.  THE  INPUT 
IS  R, H, X , S 1 , S2 .  THE  OUTPUT  I 
S  P 

120  '  PROGRAM  IS  SET  FOR  S-DECIM 
AL-OIGIT  ACCURACY  IN  P. 

125  !  ELLCV  USES  ERFC  WITH  9-DIG 
IT  RELATIVE  ACCURACY. 

130  »  SOURCES :  NHL  REPORT  *  1710 
, AUG. 1960.  MATH  OF  COMP.  OCT 
.  1961,  PP  375,382. 

135  !  “ELLCV"  CONSTRUCTED  IN  COL 
LABORATION  WITH  ALFRED  H.  MO 
RRIS  . 

140  OPTION  BASE  1 

145  DIM  P9<3>,Q9<2>,R9<5>,S9<4>, 
X<43> , Y<43> 

150  P9<1>=. 316652890658  0  P9<2>= 
1.72227577039  0  P9<3>=21.385 
33223/8 

155  Q9< 1>=7. 8437457083  §  Q9<2>=1 
8  9522572415 

1 60  R9( 1 >=4 . 31 87787405E-5  0  R9<2 
>=.56316961891  8  R9<3>=3  031 
7993362 

165  R9<4>=6. 8650184849  8  R9<5>=7 
.3738883116 

178  $9<i >=5.3542167949  8  S9<2>=1 
2  795529509 

175  S9<3>=15  18490819  8  S9<4>=7. 
3739608908 

190  X<1>=. 233619186083  8  X<2>=.6 
61209386466  8  X< 3 >=. 93246951 
4203  8  X<4>=  183434642496 

195  X<5>=. 525532409916  0  X<6>=.7 
96666477414  8  X< 7>= . 96028985 
6498  8  X<8>= . 125233408511 

200  X<9>=. 367331498998  8  X<10>=. 
587317954287  8  X< 1 1 >=. 769902 
674194  8  X<12>= .90411725637 


265  X' 13>=  .  98156063424?  8  Ji(14)  = 
9.50125098376E-2  8  X<15>=.28 
1603550779  0  X(16>= . 45801677 
7657 

210  X< 17>=. 617876244483  0  X<13>= 
.755484408355  0  Xvl9>=. 86563 
1202388  8  X<28>=. 94457582307 

3 

215  X<21>=. 989498934992  0  X<22>= 
7. 6526521 1335E-2  0  X<23>=.22 
7785351142  0  X ( 24 >=. 37370688 
8715 

220  X(25>= . 510867801951  0  X(26>= 
636053638727  0  X < 27 >=. 74633 
190646  8  X.28>=  839116971322 

225  X<29>=: 912234428251  8  X<30>= 
.953971927278  0  X(31>=. 99312 
8599185 

230  X<32>=6. 4056A928626E-2  8  XC3 
3>=. 191118867474  8  XC34>=.31 
5042679696  0  X<35>= . 43379350 
7626 

235  X<36>=. 545421471389  0  X<37>= 
648893651937  0  X<38>= . 74012 
4191579  0  X<39>=. 32000198597 

4 

240  X(48>=. 836415527004  0  X(41>= 
938274552083  0  X(42>= . 97472 
3555971  0  X<43>=. 99518721999 
7 

245  Y(l) =4679 13934573  0  Y<2>=.3 
60761573048  8  Y<3>= . 17132449 
2379  0  Y<4>=. 362683783378 

250  Y(5>= . 313706S45878  8  Y(6>=.2 
22381034453  0  Y<7>= . 10122853 
629 

255  Y<8>=. 249147845813  0  Y<9>=.2 
33492536538  0  Y< 1@>= . 2031674 
26723  0  Y(ll>=. 166078328543 

266  Y < 1 2  >  = . 106939325995  0  Y<13>  = 
4 . 71753363365E-2  0  Y(14>=  18 
9458610455  0  Y< 15>= . 18260341 
5845 

265  Y< 16>=. 169156519395  0  Y<17>= 
.149595988817  0  Y<18>=. 12462 
S971256  0  Y<19>=9. 5158511682 
5E-2  .  . 

270  Y<20>=6 . 22535239386E-2  0  Y<2 
1  >=2 .7152459411 8E-2  8  Y<22>  = 

. 152753387131  0  Y<23>=. 14917 
2985473 
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275  Y<24>=  142996109318  0  Y<25>= 
.131688638449  0  Y<26>=. 11819 
4531962  0  Y<27>=. 18193811981 
7 

280  Y< 28)  =8 . 3276741 5767E-2  0  Y<2 
9)=6 . 26720483341E-2  0  Y<30>= 
4 . 06814298004E-2  0  Y<31)=17 
6140871392E-2 

285  Y<32>=. 127938195347  0  Y<33>= 
.125837456347  0  Y<34>=. 12167 
9472928  0  Y<35>=. 11550566805 
4 

290  Y<36>=  1874442781 16  0  Y<37>= 
9 . 76 18652 104 1E”2  0  Y<38>=.08 
6190161532  0  Y<39>=7. 3346481 
41UE-2 

295  Y<40)=5  92985849154E-2  0  Y<4 
1 >=4 .42774388174E-2  0  Y<42>= 
2 . 85313886289E-2  0  Y<43>=.01 
23412298 

300  ft=4 . 892  0  91=5.387  0  82=3.87 
75  0  03=5. 16 

305  Bl=. 564139583548  0  8=1.41421 
356237  0  82=29.019769  0  !  B2 
=01401 

310  23= . 000001 4S1 4S2  0  P=0 
315  IF  R*R<=Z3  THEN  RETURN 
320  H2=H4H-*-K4K 
325  D=H0X(S1,S2) 

339  T=R-01 40 

335  !  PROCEED  TO  SEE  IF  P=0  OR  1 

340  IF  T <0  THEN  355 
345  IF  T4TCH2  THEN  355 
350  P=1  0  RETURN 

355  H8=0BS<H>  0  K8=HBS<K> 

360  IF  R— H8+04S1 <=0  THEN  RETURN 
365  IF  R-K8+04S2 <=0  THEN  RETURN 
370  S9=SQR<H2>  0  IF  S1#S2  THEN  3 
80 

375  H8=S0  0  K3=0  0  IF  R>04S1<=H8 
THEN  RETURN 

380  IF  S0<=R  THEN  435  ELSE  O=<S0 
-R)/D 

385  IF  R4R4EXP<’-<  .54040)  ><=Z3  TH 
EN  RETURN 

390  IF  S0<R+814HIN<S1,S2>  THEN  4 
35 

395  IF  H8*K8=0  THEN  435 
400  H9=H6/S1  0  K9=K8/'S2 
405  0=H94H9+K94K9  0  IF  D<=B2  THE 
N  435  ELSE  Z2=R/S2 
410  9=S2''S  1  0  Q 1  =Q4Q  0  F=Q14H94H 
9+K94K9  0  Z1=Z24Z24F/D 


415  Z=D-Z1-B2  0  IF  Z<0  THEN  425 
426  IF  Z4Z-44Z14B2>=0  THEN  RETUR 
N 

425  T 1=H84H8+Q14K84K8  0  Z8=B24S1 
4S14T1/H2  0  R2=R4R  0  Y=H2-R2 
-Z8  0  IF  Y <0  THEN  435 

439  IF  Y4Y-44R24Z8>=0  THEN  RETUR 
N 

435  Z8=0  !  FIND  LIMITS  OF  INTEGR 
0TION. 

440  Z=K8+034S2  0  H3=K8-034S2  0  S 
0=S1  0  S9=S2 

445  Z=R-Z  0  D1=0  0  IF  Z>0  THEN  D 
1=SQR<Z/R> 

459  IF  H3>  =0  THEN  455  ELSE  H5=9 
0  E3=1-D1  0  GOTO  460 

455  E3=SflR(l-H3/R>-Dl  0  H5=l 

460  IF  Z8#0  THEN  470 

465  Z8=l  0  F=E3  0  T=D1  0  Z=H8+03 
*S1  0  H6=H5  0  H3=H8-034S1  0 
GOTO  445 

470  IF  F>=E3  THEN  480  ELSE  E3=F 
475  D1=T  0  S9=S1  0  Z8=H8  0  S0=S2 
0  H8=K8  0  K8=Z8  0  H5=H6 
480  E3= . 5*E3 

485  *  G0USSI0N  INTEGR0TION  BEGIN 

S 

498  N=E3*R*<  34^S0+l/< ,025*0BS<R 
-K8)+5*S9>> 

495  Z2=R/<B*S9)  0  R8=R/'<B*S8>  0 
HQ=HQs  <B*S8>  0  K8=K8''CB*S9> 
500  IF  N<2 . 75  THEN  505  ELSE  J=31 
0  HI = 12  0  GOTO  530 
505  IF  NC1.35  THEN  510  ELSE  J=21 
0  N 1 = 1 0  0  GOTO  530 
510  IF  N < . 75  THEN  515  ELSE  J=13 
0  N1 =8  0  GOTO  530 
515  IF  NC.35  THEN  520  ELSE  J=7  0 
N 1 =6  0  GOTO  538 

520  IF  N< . 15  THEN  525  ELSE  J=3  0 
Nl=4  0  GOTO  538 
525  J=0  0  NT  =3 
530  Z=0  0  Z3=0  0  Y±81*E3*R8 
535  G3=0  0  H9= l .9999999792  0  K9= 
.0000000104  0  H=N1+N1  0  I=-N 
1 

540  IF  K8=8  THEN  635 
545  FOR  L= 1  TO  M 
550  IF  1=0  THEN  1=1 
555  T=E3*<SGN<I)tX<J+0BS< I))+ 1)+ 
0 1 

560  T9=TtT  0  T1=R8*< 1-T9)  0  T2=H 
8-T1  0  T4=EXP<-<T2ST2>> 
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565  IF  H8#0  THEN  576  ELSE  T4=T4+ 
T*  §  GOTO  589 
578  IF  H5#8  THEN  580 
575  T2=H8+T1  0  T4=T4+EXP<-< T2*T2 

>  > 

580  IF  Z#0  THEN  605 
585  Z1=Z2*T*SQR<2-T9;>  0  K1=K8-Z1 
598  IF  RBS<K1XA2  THEN  610 
595  IF  K1>0  THEN  625 
680  Z=1 

605  K5=H9  0  GOTO  615 
618  K5=FN0<K1> 

615  K1=K8+Z1  0  IF  K 1 <R2  THEN  K5= 
K5-FN0<K1 )  ELSE  K5=K5-K9 
620  G3=G3+K5*T4*T*Y<J+ABS<i;> 

625  1=1+1  0  NEXT  L 
630  P=Y*G3  0  RETURN 
635  FOR  L=1  TO  H 
640  IF  1=0  THEN  1=1 
645  T=E3*<SGN< I >*X< J+ABS< I > )+ 1 >+ 
D1 

650  T9=T*T  0  T1=R8*<1-T9)  0  T2=H 
8-T 1  0  T4=EXP<-<T2*T2)) 

655  IF  H8t0  THEN  665 
668  T4=T4+T4  0  GOTO  675 
665  IF  H5#0  THEN  675  ELSE  T2=H8+ 
T 1 

670  T4=T4+EXP<-<T2*T2>> 

675  IF  Zt0  THEN  690 
680  K1=Z2*T*SQR<2-T9>  0  IF  Kl<fl2 
THEN  K5=2*<1-FN0<K1>>  0  GOT 
0  695 
685  Z=1 
690  K5=H9 

695  G3=G3+K5*T4*T*Y<j+ABS<I>> 

700  1=1+1  0  NEXT  L 
705  P=Y*G3 
710  RETURN 

715  REM  CODY  FOR  FN0=K3=ERFC<K 1 > 
728  »  IF  K 1 <=-A2  THEN  FN0=H9 . 

725  !  IF  K1>=A2  THEN  FN0=K9 . 

730  OEF  FN0<K1  > 

735  IF  ABS<K 1 ) >  5  THEN  755 
748  K4=K1*K1 

745  FNO=l-KI*<  <P9< 1 )*K4+P9<2> ) *K 
4+P9<3>  <K4+Q9< 1 ) >*K4+Q9<2 

)) 

750  GOTO  775 
755  K4=HBS<K1> 

760  K5=<<<R9<1>*K4+R9<2>>*K4+R9< 
3>>*K4+R9<4>>*K4+R9<5> 


765  K5=K5*EXPO<Kl*Kl  >  >✓<  <  <<K4+S 
3< 1>  >*K4+S9<2>  >*K4+S9<3) >*K4 
+S9<4>  > 

770  IF  K1<0  THEN  FN0=2-K5  ELSE  F 
N0=K5 

775  FN  END 

730  X7=0  !  CRUTCHER  TfiBLE  CHECK 
785  INPUT  R,H,k,Sl,S2 

799  GOSUB  100 

795  PRINT  R;H;K;S1;S2 

800  R4=R 
805  X7=X 7+ . 1 
818  R=R4*X7 
815  GOSUB  310 

320  IMAGE  DD.DD,2X,D  DDDDDD 
825  PRINT  USING  820  ;  X7;P 
830  IF  P<=. 999999  THEN  805  ELSE 
BEEP 
835  END 
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16*  !  THIS  PROGRAM  IS  CALLED  “EL 
LCV3"  IT  SUPPLIES  THE  ELLIPT 
I CAL  COVERAGE  FUNCTION  PC R, H 
i  K  ,  S 1 S  2 ) 

165  !  P  DENOTES  THE  PROBABILITY 
OF  A  SHOT ,  NORMALLY  DISTRIBU 
TED  WITH  MEAN  C8, 8)  AND  STAN 
OARO 

116  !  DEVIATIONS  S1,S2  IN  THE  X 
Ar-tO  Y  DIRECTIONS,  RESPECTIVE 
LY, FALLING  IN  A  CIRCLE  IN  TH 
E  XY-PLANE 

115  !  OF  RADIUS  R  AND  CENTERED  A 
T  C  H ,  K  >  .  THE  OUTPUT  IS  P. 

120  !  THE  PROGRAM  IS  SET  FOR  3-D 
ECIMAL-DIGIT  ACCURACY  IN  P. 

125  !  ELLCV3  USES  ERFC  WITH  9-DI 
GIT  RELATIVE  ACCURACY. 

1 30  !  SOURCES"  NWL  REPORT  #1710, 
AUG. 1960  MATH  OF  COMP.  OCT 
1961,  PP.  375-332. 

135  i  “ELLCV3"  CONSTRUCTED  IN  CO 
LLABORAT I  ON  WITH  ALFRED  H.  M 
ORRIS 

149  OPTION  BASE  1 

145  DIM  P9C3)  ,  Q9C2),,  R9<5)  ,  S9C4)  , 
XC21),Y<21> 

150  P9C 1  )=  . 316652890658  0  P9C2)= 
1.72227577039  0  P9<3)=21.385 
3322373 

155  Q9C1 >=7.8437457033  0  Q9<2>=1 
8.9522572415 

160  R9< 1 >=4 . 31 87787405E-5  9  R9C2 
>=  56316961891  9  R9C3)=3.031 
7993362 

165  R9C4)=6. 8650184849  0.  R9<5)=7 
.3738883116 

176  S9C1 '=5. 3542167949  0  S9<2)=1 
2.795529509  0  S9<3)=25  18490 
319  0  S9<4)=7. 3739608908 

175  X(l)=. 238619186083  0  X<2)=.6 
61209336466  0  X<3)= . 93246951 
4203  0  X<4>=. 183434642496 

180  X<5)= . 525532409916  0  X<6)=.7 
96666477414  @  X<7)= . 96028985 
6498  e  X<3)=. 125233408511 

185  XO>  =.36783 1498998  0  X<10)=. 
587317954287  0  X< 1 1 )=. 769902 
674194  0  X<12>=. 90411725637 

190  XC 1 3) =. 96 1 560634247  0  XC14>= 
9 . 50125098376E-2  0  XC15)=.28 
1603550779  0  XC 16)= . 45801677 
7657 


195  XC 17 >=. 617876244403  9  X  C  1 8 > = 
755404403355  9  XC 19)=  36563 
1202388  0  XC20)=. 94457502307 
3 

200  XC 21 >=.989480934992 
205  VC 1>=. 4679 13934573  0  YC2>=.3 
60761573048  9  Y C 3) = . 1 7 1 32449 
2379  0  Y<4)=. 362683783378 
210  YC5)=. 313796645378  9  Y<6)=.2 
22381034453  8  YC7>= . 10122853 
629 

215  YC3)=. 249147945813  9  Y<9>=.2 
334S2536538  0  YC 10)= . 2031674 
26723  9  Y < 11)  =  . 160073323543 
220  Y< 12 >=. 106939325995  0  Y<13)= 
4 . 71 753363865E— 2  0  VC-14)*.  18 
9450610455  8  Y< 15)=. 18260341 
5045 

225  YC16>=. 169156519395  9  YC17>= 
149595988817  8  YC18>=. 12462 
3971256  8  YC 19)=9 . 51 5851 1682 
5E-2 

230  Y < 20) =6 . 22535239386E— 2  0  YC2 
1  )  =2 . 7152.45941 13E-2 
235  0=3.291  9  Al=3. 39895  0  A2=2 . 
898  E  A3=3 . 7 

240  Bl=. 564139583543  0  B=l. 41421 
356237  0  B2=15 .2818111  0  !  B 
2=A1*A1 

245  Z3=.801*S1*S2  9  P=0 
256  IF  R*R <=Z3  THEN  RETURN 
255  H2=H*H+K*K 
260  D=MAXCS1,S2> 

265  T=R-R1*D 

279  *  PROCEED  TO  SEE  IF  P=0  OR  1 
275  IF  T <0  THEN  285 
288  IF  T*T>=H2  THEN  P=1  9  RETURN 
285  H8=ABS(H)  0  K8=ABSCK) 

290  IF  R-H8+A*SI<=0  THEN  RETURN 
295  IF  R-K8+A*S2<=8  THEN  RETURN 
300  S0=SQR<H2>  0  IF  SltS2  THEN  3 
10 

305  H8=S0  9  K8=8  0  IF  R+«*St<=H8 
THEN  RETURN 

310  IF  S0<=R  THEN  365  ELSE  D=CS8 
-R)/D 

315  IF  R*R*EXPC-C .5*0*0 >)<=Z3  TH 
EN  RETURN 

320  IF  S0<R+A1*MIN<S1,S2)  THEN  3 
65 

325  IF  H8*K8=0  THEN  365 
338  H9=H8/S1  0  K9=K8/S2 
335  D=H9*H9-HC9*K9  0  IF  D<=B2  THE 
N  365  ELSE  Z2=R/S2 
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346  Q=S2/S1  0  Q1=Q*Q  6  F=Q1*H9*H 
9+K9XK9  0  Zl=Z2XZ2XF/0 
345  Z=D-21— 62  0  IF  Z<8  THEN  355 
350  IF  Z*Z-4*21*B2>=9  THEN  RETUR 
N 

355  Tl=H8*H8+Ql*k8*K8  §  Z8=B2*S1 
XS1XT1/H2  0  R2=R*R  0  Y=H2-R2 
-23  0  IF  Y <8  THEN  365 
368  IF  Y*Y-4*R2*Z8>=0  THEN  RETUR 
N 

365  28=8  !  FIND  LIMITS  OF  INTEGR 
ATION 

378  Z=K8+A3*S2  0  H3=K8-A3*S2  0  S 
0=S1  0  S9=S2 

375  2=R-2  0  01=8  0  IF  2>0  THEN  D 
1 =SQR  CZ/R  > 

339  IF  H3>  =0  THEN  385  ELSE  E3=l- 
D1  0  H5=0  0  GOTO  398 
385  E3=SQR<1-H3/R>-01  0  H5=l 
390  IF  28#0  THEN  490 
395  28=1  0  F=E3  0  T=01  0  Z=H8+A3 
*S1  0  H€=H5  0  H3=H8-fl3XSl  0 
GOTO  375 

406  IF  F>=E3  THEN  410  ELSE  E3=F 
485  D1=T  0  28=H8  0  S0=S2  0  S9=S1 
0  H8=K8  0  K8=Z8  0  H5=H6 
410  E3= . 5XE3 

415  »  GAUSSIAN  INTEGRATION  BEGIN 

S 

429  N=E3*Rt< . 34/S0+l/< . 825XABSCR 
-K8  >+5XS9  >  !> 

425  Z2=R/<B*S9>  0  R8=R/<BXS0>  0 
H8=H8/<B*S0>  0  K8=K8/<BXS9> 

430  IF  N<2  THEN  435  ELSE  J=13  0 
Nl=8  0  GOTO  458 

435  IF  N< . 675  THEN  448  ELSE  J=7 
0  Ni=6  0  GOTO  450 
440  IF  N< . 5  THEN  445  ELSE  J=3  0 
N 1 =4  0  GOTO  450 
445  J=0  0  N1 =3 
450  7=0  0  Z3=0  0  Y=B 1 XE3XR8 
455  K9=. 000015  0  H9=l. 99997 
460  G3=0  0  M=N1+N1  0  I=-N1 
465  IF  K8=8  THEN  560 
470  FOR  L= 1  TO  M 
475  IF  1=8  THEN  1=1 
488  T=E3*<SGN<I>*X<J+ABS<I>)+1>+ 
01 

485  T9=T*T  0  71=R8*<1-T9>  0  T2=H 
8-T1  0  T4=EXP<-<T2*T2>  > 

490  IF  H8#0  THEN  495  ELSE  T4=T4+ 
T4  0  GOTO  585 
495  IF  H5#9  THEN  595 


586  T2=H8+T 1  0  T4=T4+EXP<-<T2*1"2 

>> 

505  IF  2#0  THEN  530 

518  Z1=Z2*T*SQR<2-T9>  0  K1=K8-Z1 

515  IF  K1CA2  THEN  535 

528  IF  Kl>8  THEN  550 

525  2=1 

538  K5=H9  0  GOTO  540 
535  K5=FNA<K1> 

540  K1=K8+Z1  0  IF  K1'!A2  THEN  K5= 
K5-FNA<K1 >  ELSE  K5=K5-K9 
545  G3=G3+K5*T4*T* Y  <  J+ABS  < I >  > 

558  1=1+1  0  NEXT  L 
555  P=Y*G3  0  RETURN 
560  FOR  L=1  TO  M 
565  IF  1=8  THEN  1=1 
578  T=E3*<SGN< I >*X< J+ABS< I > >+l >  + 
01 

575  T9=T*T  0  T1=R8*<1-T9>  0  T2=H 
8-T 1  0  T4=EXP<-<T2*T2>> 

580  IF  H8#0  THEN  590 
585  T4=T4+T4  0  GOTO  608 
590  IF  H5#8  THEN  600  ELSE  T2=H8+ 
T 1 

595  T4=T4+EXP<-<T2*T2>> 

660  IF  Z#9  THEN  615 
605  K1=Z2*T*SQR<2-T9>  0  IF  K1<A2 
THEN  K5=2X< 1 -FNA<K1 > )  0  GOT 
0  628 
618  2=1 
615  K5=H9 

620  G3=G3+K5*T4*TXY< J+ABSC I > > 

625  1=1+1  0  NEXT  L 
638  P=YXG3  0  RETURN 
635  REM  CODY  FOR  FNA=K3=ERFC  <  K 1 ) 
-9-DIGITS. 

640  !  IF  K1<=-A2THEN  FNA=H9 .  IF 
K1  >  A2  THEN  FNA=K9 
645  DEF  FNA<K1> 

650  K5=0  0  IF  ABS  < K 1 > > . 5  THEN  67 
0 

655  K4=K1*K1 

660  FNA=1-K1X<  <P9< 1 )XK4+P9<2>  >XK 
4+P9<3>  >/<,  <K4+Q9<  1  >  >*K4+Q9<2 

)) 

665  GOTO  690 
670  K4=ABS<K1> 

675  K5=<  <  <R9< 1 >Xk4+R9(2>  >XK4+R9( 
3) )XK4+R9<4>  >XK4+R9<5) 

680  K5=K5/<  <  <  <K4+S9< 1 >  >*K4+S9C2> 
>*K4+S9<3>>*K4+S9<4>>*EXP<-< 
K1XK1 >  > 

685  IF  K 1 <0  THEN  FNA=2-K5  ELSE  F 
NA-k5 

698  FN  END 
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POLYCV 

Input:  (Data  Statement  for  HP-9845)  P8,  P9,  Mx,  My>C,  Sx.  Sy,  N 

(Data  Statement  for  HP-85)  P8,  P9,  Ml,  M2,C,SI,S2,N 

(Mx,  My)  or  (Ml,  M2)  =  mean  of  the  normal  distribution 
c  =  correlation  coefficient  (contained  in  C) 

Sx,  Sy  or  SI,  S2  =  standard  deviations  ax  and  oy 

N  =  1,  K  =  3  =*  probability  desired  over  a  single  angular  region  Al.  Vertex  of  A1  is 
always  given  by  (Xj,  yj).  Points  (x2,  y2),  (x^,  y3)  are  given  in  counterclockwise 
order  about  the  vertex.  POLYCV  sets  K  =  3. 

N  =  K  >  3  ■*  probability  desired  over  a  polygon  E.  Vertices  are  specified  in  counter¬ 
clockwise  order.  POLYCV  sets  K  =  N. 

P8  is  used  to  specify  how  the  coordinates  of  the  points  defining  E  or  Al  are  stored.  If 
P8  =  0,  then  Xj,  yb  —  xK,  yK  are  stored  consecutively  in  data  statements  im¬ 
mediately  following  the  initial  data  statement  above.  POLYCV  then  stores  Xj  in 
array  element  X(J)  and  yj  in  array  element  Y(J).  If  P8  =£  0,  then  it  is  assumed  by 
POLYCV  that  the  points  are  already  stored  in  arrays  X(*),  Y(*).  This  is  useful  if 
the  points  are  machine  generated. 

P9  is  used  to  specify  what  output  is  desired  as  indicated  below. 

Output:  P,  A,  W,  1 1  are  always  given  as  part  of  the  output  if  N  >  3.  For  N  =  1 ,  P  ar.d  1 1  are 

given. 

P  =  probability  from  a  normal  correlated  bivariite  distribution  over  an  arbitrary  polygon 
E  or  an  angular  region  A 1 . 

A  =  area  of  E.  If  N  =  1 ,  then  A  is  set  to  zero. 

W  =  winding  number  of  E.  For  simple  polygon s,  W  *  I .  If  N  -  1,  W  is  set  to  zero.  W  is 
contained  in  W1 . 

II  =  Error  and  Information  Parameter: 

II  =0,  output  acceptable. 

II  55  -1,  angular  region  Al,  N  =  I,  may  not  be  well-defined.  The  angular  measure 
of  A 1  is  close  to  0  or  2».  A  result  for  P  is  given. 

11  =  1,  angular  region  Al,  N  =  1,  is  not  well-defined,  i.e.,  at  least  one  of  the  two 
points  (x2,  y2),  (x3,  y3)  is  too  close  to  (Xj,  y,).  Input  unacceptable.  P  is  set 
to  -5. 
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11  =  2,  two  consecutive  segments  of  E  overlap.  Output  is  O.K. 

II  =3,  the  correlation  coefficient,  c,  contained  in  C,  does  not  satisfy  |c| 
Input  is  unacceptable.  P  is  set  to  -5. 

P9  =  0,  no  additional  output  besides  the  above  is  given. 

P9  =  1,  the  x,  y  coordinates  of  the  vertices  of  E  or  the  points  of  A 1  arc  listed. 

P9  <  0,  a  plot  of  E  or  A 1  is  given. 

P9  >  1 ,  both  a  listing  of  the  x,  y  coordinates  and  a  plot  of  E  or  A1  are  given. 
Accuracy:  P  given  correctly  to  approximately  9-decimal  digits. 
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P8-  0 

P9-  2  Mx=- . 5 

My  = 

6  C 

=  0  Sx 

= 

X< 

1 

>-  0 

Y  < 

1  >- 

6 

X< 

? 

>«  1 

Y  < 

2  )=■ 

1 

X< 

3 

>«  8 

Y  < 

3  >- 

'"t 

X< 

4 

>  — 1 

Y  < 

4  )  - 

2 

X< 

5 

>-  0 

Y< 

5  >- 

1 

X< 

6 

>  — 2 

V  < 

6  >- 

0 

XC 

7 

>«- 1 

Y  < 

7  >« 

0 

X( 

8 

>»  0 

Y  < 

8  >  — 

1 

X< 

9 

>■-1 

Y  < 

9  >  — 

1 

X< 

10 

>-  0 

Y  < 

10  >■ 

-2 

x< 

11 

>«  1 

Y  < 

11  >■ 

0 

p* 

.26421465361  ft- 

4.5 

W- 

1  11- 

0 

11 


P-. 26421465361 


fl-4 
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P8*  0 

P9*  2  fix* 

.3  My*  1  C 

*  0  Sx*  1 

X< 

1 

>*- 3 

Y  < 

l  )■  4 

X< 

2 

>*-3 

Y  < 

2  )  — 3 

X< 

3 

)  *  4 

Y  < 

3  >*-3 

X< 

4 

>*-l 

Y  < 

4  >*  2 

XC 

3 

)*  2 

Y  < 

3  )  — 1 

X< 

6 

)*  4 

Y< 

6  >»  1 

x< 

7 

)«  0 

Y  < 

7  >*  1 

x< 

8 

)*  0 

Y  < 

8  )—  2 

x< 

9 

>*- 2 

Y  < 

9  )*  0 

x< 

10 

)»  1 

Y  < 

10  )■  0 

p- 

.40616638044 

0*  27  W*  2 

11*  2 

P». 40616638044 


Sy=  3  N. 


R 
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P8 

*  1 

P9*  2  Mx*  0  My*  0 

C* 

0  Sx*  3  Sy*  2  N*  12 

X< 

1 

>*  7 

Y  < 

1  >-  3 

X< 

2 

>*  6.33012701885 

YC 

2  >*  5.5 

XC 

3 

>*  4.49999999987 

Y  < 

3  >■  7.33012701885 

X< 

4 

>*  1.99999999988 

Y  < 

4  >=  7.99999999981 

X< 

5. 

>*-.50000000001 

Y  < 

5  >*  7.33012781861 

X< 

6 

>*-2.33012701873 

Y  < 

6  >-  5.49999999966 

X< 

7 

>*-2.9999999996 

Y< 

7  >-  2.99999999976 

X< 

8 

>*-2.33012701838 

Y< 

8  >*  .49999999999 

X< 

9 

>*-.49999999946 

Y< 

9  >*-1.33012701861 

x< 

18 

>«  2.00000800035 

Y< 

10  >*-1,99999999939 

x< 

11 

>*4.5 

Y< 

11  >*-1.33012701815 

x< 

12 

)*  6.3301270185 

Y< 

12  >*  .50000000074 

P*  .56768766703  A*  74.9999999898  W*  1  II*  8 


P». 56768766703  fl-74 . 9999999898  H-l 

Y 


59 


"■*./ «  •**  v 


*  «%  •'»  *«.•  <“•  '**V  »• 


-  •  >  *  •  •  T  »  *  7  m  -  *  >  * 

'  4.*  *  .  •„*  «L  *1  "  v  M  ^  «, 


-  .f  *  •  -  .*  « 

'•  Ofv'  v*  VI  ■/ 
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P8*  0  P9*  2  Mx*  0  My*  0  C*  .  3  Sx  =  1  Sy*  2  N*  i 


X<  1  >*  1 
X<  2  >*  1.5 
X<  3  >*  2 


Y(  1  )*  1 
Y<  2  )  *  2 
Y <  3  >*  .5 


P*  .95617907488  fl*  0  W*  0  II*  0 


P-. 95617907488 


R=0 
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INPUT =  0  2  .5  1  0  1  3 


X<  1  >»-3 
X(  2  >=-3 
X<  3  >=  4 
X<  4  )=-l 
X<  5  >*  2 
X<  6  >=  4 
X<  7  >=  0 
X<  8  >=  0 
X<  9  >=-2 
X<  10  )=  1 


Y<  1  >=  4 

«  §  £3 

Y<  4  )=  2 
Y <  5  >=-l 
Y<  6  >=  1 
Y<  7  >  =  1 
Y  <  8  >=-2 
Y<  9  >-  0 
Y<  10  )=•  0 


10 


.40616638001  27  2  2 
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INPUT:  1  2  0  0  0  3  2  1 

X  <  1  >  =  7  Y  <  1  >  =  3 

X(  2  >  =  6.33012701833 

Y <  2  )=  5.5 
X<  3  >  =  4.50080000001 

Y  <  3  >  = 

7 . 33012701893 
y,<  4  >  =  2  Y<  4  )  =  8 

X(  5  >  =  - .  5  Y<  5  >  = 

7  33012701893 
X<  6  >=-2.33012701893 

Y  <  6  >  = 

5  50000000001 
X<  7  >=-3  Y<  7  >=  3 

X<  8  >=-2.33012701893 

Y <  8  >=  .5 

X<  9  >=-.50000000081 

Y  <  9  >  = 

-1  33012701893 

X<  10  >=  2  Y <  10  >=-2 

X<  11  >=  4.5  Y<  11  >= 

-1 .33012701893 
X<  12  >=  6.33012701893 
Y<  12  >= 

.49999999999 


.567687666698  75  1  0 


NSWCTR  83-13 
POLYCV-HP  9845 

109  !  THIS  PROGRAM  IS  CALLEfl  "POLYCV".  IT  SUPPLIES 

THE  PROBABILITY  P  OF  A  SINGLE  SHOT,  NORMALLY 
105  !  DISTRIBUTED  IN  THE  XY-PLANE  WITH  MEAN  <Mx,My> 

STANDARD  DEVIATIONS  Sx.Sy  AND  CORRELATION 

110  !  COEFFICIENT  C,  FALLING  IN  AN  ARBITRARY  POLYGON 

E  OR  IN  A  SEMI-INFINITE  ANGULAR  REGION  A1  IN  THE 
115  ■  PLANE.  E  IS  SPECIFIED  BY  K  POINTS  <X(J>,Y<J>>.  FOR  J  = 

1  TO  .  K .  A1  IS  SPECIFIED  BY  3  POINTS  <X(J),Y<J>>  WITH 
120  !  N= 1  AND  THE  VERTEX  OF  A1  GIVEN  BY  <X(1>,Y<1>>.  THE 

POINTS  MUST  BE  GIVEN  IN  COUNTERCLOCKWISE  ORDER. 

125  !  THE  INITIAL  INPUT  IS: P8, P9, Mx, My, C, Sx, Sy, N.  IT  IS  STORED 

AS  DATA  BEGINNING  AT  2000.  IF  N*K>=3 
130  !  THEN  THE  INPUT  SPECIFIES  A  POLYGON  E.  IF  N=1  (WITH 

K=3) ,  THEN  AN  ANGULAR  REGION  A1  13  GIVEN. 

135  !  IF  P8*0,  THEN  THE  COORDINATES  X(J),Y(J>,  FOR  J=1  TO  K, 

ARE  STORED  IN  DATA  STATEMENTS  IMMEDIATELY  FOLLOWING 
140  !  THE  INITIAL  INPUT  DATA  STATEMENT.  POLYCV  STORES  THEM  IN 

ARRAYS  X<*>, Y<*>.  IF  P8§0,  THEN  IT  IS  ASSUMED  THE  VERTEX 
145  !  COORDINATES  ARE  ALREADY  STORED  IN  THE  ARRAYS  X(*),Y<*>. 

THIS  IS  USEFUL  IF  THE  VERTEX  COORDINATES  ARE  MACHINE 
150  !  GENERATED.  P  IS  GIVEN  TO  NINE  DECIMAL-DIGIT-ACCURACY. 

155  !  IF  P9  >*  1  THEN  LIST  X(J),Y<J>.  IF  P9>1  OR  <0  THEN  PLOT  E 

OR  fll.  IF  P9*0  THEN  NO  PLOT  OR  LIST. 

160  !  THE  OUTPUT  IS:  P,A,W,I1,  WHERE  A  CONTAINS  THE  AREA 

OF  E,  W1  CONTAINS  THE  WINDING  NUMBER  W  OF  E,  AND  11 
165  !,  IS  AN  ERROR  INDICATOR.  IF  11*0  OR  2  THEN  THE  OUTPUT  IS 

ACCEPTABLE.  IF  II  *1  OR  -1,  THE  ANGULAR  REGION  A1  MAY  NOT 
170  !  BE  WELL-DEFINED.  IF  11*1  THE  VERTEX  AND  ONE  OF  THE  OTHER  TWO 

POINTS  MAY  BE  TOO  CLOSE  TO  EACH  OTHER.  IF  II  — 1,  THEN  THE 
175  '  MEASURE  OF  A1  IS  CLOSE  TO  0  OR  2PI. 

130  !  IF  11*2  THEN  TWO  CONSECUTIVE  SIDES  OF  E  OVERLAP.  OUTPUT  IS  O.K 

IF  11*3  THEN  C*C>*1.  IF  11*1  OR  3  THEN  P  SET  TO  -5. 

185  !  SOURCES:  NSWCxDL  REPORT  #3886 , SEPT . 1 978 .  NSWC/DL  REPORT#80- 1 66 

JUNE, 1980.  SIAM  JN. SCI. STAT. COMPUT. , JUNE  1980, PP. 179, 186. 

190  !  SIAM  JN. SCI. STAT. COMPUT. , DEC, 1982,  PP.  ?. 

195  !  X(J)  AND  Y(J>  DIMENSIONED  AT  90.  IF  MORE  POINTS  ARE 

NEEDED  TO  SPECIFY  E  THEN  MAKE  CHANGES  AT  LINE  200. 

200  DIM  X(90>,Y(90),Ul<15> 

205  U1 ( 1 >*4.08335517232E-7 
210  U1  <2>— 9.  7186486416E-6 
215  U 1 < 3 >* 1 . 057875 7448 1 E-4 

220  U 1 ( 4 ) — 7. 04260243309E-4 
225  U1 <5>*3. 24944543171E-3 
230  U1 <6>  — 1 . 12323532148E-2 
235  U1 (7>*3. 09199295521E-2 
240  U1  <8>  — 7.  149098378E-2 

245  IJ1  (9)-.  145060043403 
250  U1  ( 10)  — .  265638206366 

255  U1 ( 1 1 >*.442851899329 
260  Ul<12>  — .  666626670511 

265  U1 < 13)*. 886223733187 
270  Ul < 14)*-. 999999999776 
275  U1(1S>*. 886226924931 
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280 
285 
290 
295 
300 
305 
310 
315 
320 
325 
330 
335 
340 
345 
350 
355 
360 
365 
370 
375 
380 
385 
390 
395 
400 
405 
410 
415 
420 
425 
430 
435 
440 
445 
450 
455 
460 
465 
470 
475 
480 
485 
490 
495 
500 
505 
510 
515 
520 
525 
530 
15 


29=1.12837916709 
81=2. 71E-19 
82=1 . 34E-6 


84  =  7. 31  IE-4 
82=19.201924 

Rl*. 564189583546  !  1/SQR<PI> 
R3=4. 9E-27 


T7=3. 14159265359 
R4*. 159154943092 
T9* 1 . 14472988585 
R6=. 28209479177 
T8*7E-12 


87=3. 2625E-1 1 

!  INPUT:  P8,P9,Mx.My,C,Sx,Sy,N,X< J>, Y< J>  <J=1  TO  K> 

RE8D  P8,P9,Mx,My,C,Sx,Sy,N 
EXIT  8LPH8 
EXIT  GR8PHICS 

print  "P3=";P8; "P9*“;P9; "Hx=":mx; "My**; My; HC=";C; "Sx=*;Sx; 
PRINT 


C7=1-C*C 

IF  C7>0  THEN  405 

11-3 

P—5 

print  -ii»-;3; -p»“;-5 

RETURN 

C7=SQR<C7> 

K-3 

IF  NOl  THEN  K-N 
FOR  J-l  TO  K 

IF  P8-0  THEN  RE8D  X<J),Y<J) 

if  P9 >= l  then  print  “x j ;*>-"; x < j >; T8B< 25) ; "Y< B; j Y c , 
NEXT  J 

IF  < P9 >  1 )  OR  <P9<0>  THEN  1560 
FOR  J-l  TO  K 

Y< J>-<Y< J)-My)'Sy 
X<J)-<<XCJ)-Mx)>Sx-C*Y<J>>-'C7 


NEXT  J 

X<K*1 )*X<  1  >  , 

Y(K*1 >»Y<  1  > 

P-0 
11-0 
8-0 
Kl-0 
K- 1 

IF  NOl  THEN  550 
X1-X<1) 

Y 1 -Y<  1  ) 

X<3>-X<l>*X<l>-X<3> 
Y<3>-Y<i>*Y<l>-Y<3> 
U-X<2)-X< 1 > 
V-Y<2)-Y< 1 > 

W-X< 1 >-X<  3> 

2-Y< l >-Y<3> 


66 


C 
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545  GOTO  590 

550  V1=Y<N+1) 

555  X 1 =X  <  N+ 1 > 

560  U=X  <  2  > -X  <  1  ) 

565  V=Y(2)-Y<1> 

570  X1=X<1> 

575  Y1*Y<1> 

580  W=X  < 1 ) -X ( N ) 

585  Z=Y  < 1 ) -Y  <  N ) 

590  D 1 =W*W+Z*Z 

5S5  IF  D1 >R3  THEN  620 

600  IF  N=1  THEN  1370 

605  N=N-1 

610  . IF  N*2  THEN  1175 

615  GuTO  580 

620  D2=U*U+V*V 

625  IF  D2>R3  THEN  665 

630  IF  N* 1  THEM  1370 

635  K*K  + 1 

640  U=X<K+1)-X1 

645  V*Y  <K+ 1  )  -Y 1 

650  D2*U*U+V*V 

655  IF  D2<  =R3  THEN  635 

660  IF  K*N- 1  THEN  1175 

665  fl=X:*<Y<K+l)-Y<N) > 

670  B1*SQR<2*D1> 

675  B2=SQR<2*D2> 

680  P2*V*W-U*Z 

685  C1=U*W+V*Z 

690  GOSUB  1390 

695  K 1  *K  1 >05 

700  L*0 

705  B=».5*<X<K)*X<'K)+Y<K)*Y<K)> 

710  IF  B >fl 1  THEN  730 

715  C2=0 

720  P 1 *05*R4-C2 

725  GOTO  1145 

730  G1*<M*X<K)+Z*Y<K) )/Bl 

735  G2*<U*X(K)+V*Y<K) )/B2 

740  H1'»<Z*X<K)-W*Y(K),  )^B1 

745  H2*<V*XOO-U*Y<K>  >  .'B2 

750  IF  PBS<P2>>2*Bl*B2*fl7  THEN  830 

755  IF  Cl<0  THEN  780 

760  IF  fiBS< P5  >  <  *T8  THEN  770 

765  IF  Gl<0  THEN  830 

770  P1*0 

775  GOTO  1145 

780  11*2 

785  IF  P2<  0  THEN  810 

790  H»H2 

795  GOSUB  1435 

800.  PI*. 5*E 1 

805  GOTO  11.45 
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810  H=H1 

815  GOSUB  1435 

820  P 1  *- < . 5*£ 1 ) 

825  GOTO  1145 
838  IF  B< =A2  THEN  1050 
835  IF  Gl<  0  THEN  900 
840  IF  G2>*0  THEN  1090 
845  G2*-G2 

850  H2=-H2 

855  IF  fiBS(H2)<=fl4  THEN  880 

860  H=-H2 

865  GOSUB  1435 

870  L». 5*E1 

875  GOTO  1060 

880  L*.5+Rl*H2 

885  GOTO  1060 

890  L=.5-R1*H1 

895  GOTO  1060 

900  Gl—Gl 

905  H 1 *-H 1 

910  IF  G2< 0  THEN  940 

915  IF  RBSCHl  ><*A4  THEN  890 

920  H»H1 

925  GOSUB  1435 

930  L*. 5*E1 

935  GOTO  1060 

940  G2—G2 

945  H2— H2 

950  IF  flBS<Hl)OA4  THEN  1815 

955  IF  ABS<H2>Ofl4  THEN  995 

960  H-Hl 

965  GOSUB  1435 

970  L».5*E1 

975  H*H2 

980  GOSUB  1435 

985  L*L- . 5*E 1 

990  GOTO  1090 

995  H»H1 

1600  GOSUB  1435 

1005  L«R1*H2-.5*<1-E1> 

1010  GOTO  1090 

1015  IF  RBS<H2X*A4  THEN  1040 

1020  H-H2 

1025  GOSUB  1435 

1030  L«.5#<1-E1>-R1#H1 

1035  GOTO  1090 

1040  L*R1*<N2-H1 > 

1045  GOTO  1090 

1050  C2»R6*<H2-H1 )-R4*<G2*H2-Gl*Hl ) 

1055  GOTO  720 
1060  P2— P2 

1065  IF  P2O0  THEN  1085 
1070  L-L-l 
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1075  R5*T7+R5 

1080  GOTO  1080 

1085  R5*R5-T7 

1090  IF  B  >*R2  THEN  1140 

1095  C3»R5 

1100  C4».5*R5 

1105  M*15 

1110  F>0 

1115  R6-H2-H1 

1120  C5-R6 

1125  GOTO  1310 

1130  P1»L+EXP<-<B+T9))«<C4-S2) 

1135  GOTO  1145 
1140  Pl^L 

1145  IF  KON  THEN  1235 
1150  IF  NOl  THEN  1190 
1155  P«RBS<P1> 

1160  IF  K 1 >0  THEN  P-flBS<l-P> 

1165  C4-W1-0 

1170  IF  RBS<K1  )<  IE-1 1  THEN  II  — 1 
1175  PRINT 

1180  PRINT  “P--;P; -Ii»-; n 

1185  GOTO  1485 

1190  P-P-Pl 

1195  K1«K1*R4 

1200  fl-.5*fl 

1205  IF  K 1 < 0  THEN  1220 

1210  H1»INT<K1+. 1> 

1215  GOTO  1225 
1220  14 1  ■  I  NT  <  K1  ♦ .  9  ) 

1225  P-P+Wl 
1230  GOTO  1470 
1235  U-U 
1240  Z-V 
1245  B1-B2 

1250  X1*X<K*1 > 

1255  Yi»Y<KM> 

1260  Y2*Y  <  K) 

1265  K“K  + 1 

1270  U»X<K>1 >-Xl 

1275  V«Y<I01  )-Yl 

1280  D2«U*U+V*V 

1285  IF  92<-R3  THEN  1265 

1290  B2«SQR<2*D2> 

1295  P-P-Pl 

1300  fl-R+Xl*(Y<K+l )-Y2) 

1305  GOTO  683 
1310  S2»U1<M>*R6 

1315  M»M-1 

1320  H2«H2*G2 

1325  H 1 *H 1 *G 1 
1330  T-H2-H1 

1335  F«F  +  B 
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1340  C6-<F#C3+T)^<16-M) 

1345  S2»S2+U1<M)*C6 

1350  IF  M«1  THEN  1130 

1355  C3-C5 

1360  C5-C6 

1365  GOTO  1315 

1370  P—5 

1375  11*1 

1380  print  •n»";i;"P*";-5  ■ 

1385  RETURN 

1390  IF  ABS<C1X»ABS<P2>  THEN  1415  ! A5=ARCTANCENT<P2/C1 > ,  <-PI,PIJ. 

1395  IF  <C1<0)  AND  <P2»0>  THEN  1425 
1400  A5«ATN<P2/C1> 

1405  IF  C1<0  THEN  A5-A5+SGN < P2 > *T7 
1410  RETURN 

1415  A5«SGN<P2>»<  .  5#T7-ATN<C1/'ABS<P2>  >  > 

1420  RETURN 
1425  A5-T7 

1430  RETURN 

1435  El-0  !  £1-ERFC<H> 

1440  C4»ABS<H) 

1445  IF  C4>*4. 4  THEN  1460 

1450  El-<<<<<<<<<<Ul<l>*C4+Ul<2>>*C4+Ul<3>>*C4«-Ul<4>>#C4+i:is5>>*C4+Ul 

<6))*C4+U1<7))*C4+U1<8))*C4+U1<9) )»C4+U1 <10>>*C4+U1<11>)*C4+U1C12) 

1455  E1-<<<E1*C4+U1<13>>*C4+U1<14) >*C4»U1 < 15) >*Z9*EXP(-H»H> 

1460  IF  H< 0  THEN  E1-2-E1 
1465  RETURN 
1470  C4-Sx*Sy*C7#A 
1475  PRINT 

1480  print  "P-“;p; "A-*;C4; "M-“;ui; "II-"; II 
1485  IF  <  P9> 1 )  OR  <P9<0>  THEN  1905 
1490  BEEP 
1495  RETURN 

1500  PLOTTER  IS  "GRAPHICS* 

1505  GRAPHICS 
1510  U»X<1> 

1515  V-U 
1520  H-YU) 

1525  Z-U 

1530  FOR  1*1  TO  K 

1535  IF  X<I)>U  THEN  U-X<I) 

1540  IF  X<IXV  THEN  V-X<I> 

1545  IF  Y< I ) >M  THEN  H-Y<I> 

1550  IF  Y< I XZ  T  EN  Z-Y<I> 

1555  NEXT  I 

1560  GCLEAR 
1565  Pl-U-V 

1570  IF  P1O0  THEN  1590 
1575  Pl-i 

1580  V-V-Pl 
1585  U-U+Pl 
1590  Yi -LGT<P1 > 

1595  X1-INT<Y1> 
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1600  V 1 “FRACT  <  Y 1  ) 

1605  Y1*10''Y1 

1610  IF  Xt<0  THEN  Y 1  *  1 0*Y 1 

1615  Xl^lO^Xl 

1620  IF  Y1 >=2  THEN  1635 

1625  T 1  * . 1 *X 1 

1630  GOTO  1670 

1635  IF  Y 1 >=4  THEN  1650 

1640  T 1  * . 2*X  1 

1645  GOTO  1670 

1650  IF  Y 1 >*5  THEN  1665 

1655  T 1 = . 4*X 1 

1660  GOTO  167fa 

1665  T1 *. 5*X1 

1670  T3=INT<V>T1)*T1 

1675  T  4=11 

1680  IF  FRACT«;T4/T1><>0  THEN  T4* < I  NT < T4/T 1 >♦ 1 > *T1 
1685  P2*U-Z 

1690  IF  P2O0  THEN  1710 
1695  P2= 1 

1700  Z=2-P2 

1705  W=W+P2 
1710  Y 1 =LGT  <  P2 ) 

1715  X1=INT(Y1) 

1720  Y 1 “FRflCT  <  Y 1  ) 

1725  Y1*10-Y1 

1730  IF  Xi<0  THEN  Y1»10*Y1 

1735  X1-10-X1 

1740  IF  Yl>=2  THEN  1755 

1745  T2=.1*X1 

1750  GOTO  1790 

1755  IF  Yl>«4  THEN  1770 

1760  T2=.2*X1 

1765  GOTO  1790 

1770  IF  Y 1 >*5  THEN  1785 

1775  T2=.4*X1 

1780  GOTO  1790 

1785  T2-.5#X1 

1790  T5»INT<Z/T2)*T2 

1795  T6-0+W 

1800  IF  FRfiCT  <T6^T2><>0  THEN  T6» < I  NT < T6/T29  +  1 > #T2 

1805  SCALE  T3-.35*P1,T4+.25*P1,T5-.25*P2,T6+.2*P2 

1810  CLIP  T3,T4,T5,T6 

1815  LAXES  T1,T2,T3,T5,5,5,4 

1820  LflXES  Tl,T2,T4,T6,3,5,2 

1825  IF  N-l  THEN  1945 

1830  MOVE  X<1>,Y<1> 

1835  FOR  I-l  '.0  N 
1840  DRAW  X< I > , Y< I > 

1845  LABEL  VAL*<I> 

1850  MOVE  X<I),Y<1) 

1855  NEXT  I 

1860  DRAM  X<1>,Y<1> 
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1863  HOVE  U*T1,Z-.02*P1 
1870  LABEL  "X“ 

1873  HOVE  T3-.02#P1,T6+.05*P2 

1880  LABEL  "V- 

1883  C8»T3 

1890  C9-T6-t.l»P2 

1893  D3«P1 

1900  GOTO  443 

1903  HOVE  C8,C9 

1910  LABEL  "P»“liVAL«<P) 

1913  HOVE.  C8*.5*D3,C9 
1920  LABEL  "A-“tVAL*<C4) 

1925  HOVE  C8+D3.C9 
1930  LABEL  “H-- tVAL»<Ul ) 

1933  DUHP  GRAPHICS 

1940  RETURN 

1943  HOVE  X<3>,Y<3> 

1930  LABEL  VAL*<3) 

1935  HOVE  X<3>,Y<3> 

1960  DRAM  X<l>,v<l> 

1963  LABEL  VAL»<1) 

1970  HOVE  X<1),YU) 

1973  DRAM  Xt2)fY<2) 

1980  LABEL  VAL«<2) 

1983  GOTO  1865 
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180  !  THIS  PROGRAM  IS  CALLED  "PO 
LYCV" .  IT  SUPPLIES  THE  PROBA 
BILITY  P  OF  A  SINGLE  SHOT,  N 
ORMALLY 

105  !  DISTRIBUTED  IN  THE  XY-PLAN 
E  WITH  MEAN  <M1 , M2 >, STANDARD 
DEVIATIONS  S1,S2  IN  THE  X  A 
NO  Y 

119  !  DIRECTIONS, RESPECTIVELY, 
CORRELATION  COEFFICIENT  C, 

115  »  FALLING  IN  AN  ARBITRARY  PO 
LYGON  E  OR  A  SEMI-INFINITE  A 
NGULAR  REGION  A1  IN 

120  !  THE  XY-PLANE.  E  IS  SPECIFI 
ED  BY  K  POINTS  <X< J) , Y< J) ) , J 
=1  TO  K.  fll  IS  SPECIFIED  BY 
3  POINTS 

125  !  <X< J) , Y< J) )  WITH  THE  VERTE 
X  OF  A1  GIVEN  BY  <X<1>,Y<1>> 
AND  THE  POINTS  GIVEN  IN  COU 
NTER 

130  !  CLOCKWISE  ORDER.  THE  INITI 
AL  INPUT  IS : P8, P9, Ml , M2, C, SI 

,S2,N. 

135  !  IT  IS  STORED  AS  DATA  BEGIN 
NING  AT  1085. 

140  !  IF  N=K>=3  THEN  THE  INPUT  S 
PECIFIES  E.  IF  N=1 (WITH  K=3> 

, THEN  A1  IS  GIVEN. 

145  !  IF  P8=0,  THEN  THE  COORDINA 
TES  X<J),Y<J),  FOR  J=1  TO  K, 
ARE  STORED  IN  DATA  STATEMEN 
TS  IM- 

158  !  MEDIATELY  FOLLOWING  THE  IN 
ITIAL  DATA  STATEMENT.  POLYCV 
STORES  THEM  IN  ARRAYS  X<*>, 

Y<*>. 

155  1  IF  P8#0,  THEN  IT  IS  ASSUME 
D  THE  VERTEX  COORDINATES  ARE 
ALREADY  STORED  IN  ARRAYS  X< 

*>,Y<*>. 

160  !  THIS  IS  USEFUL  IF  THE  VERT 
EX  COORDINATES  ARE  MACHINE  G 
ENERATED.  IF  P9>  =  1  THEN  LIST 

X<J),Y<J>. 

165  '  IF  P9> 1  OR  <0  THEN  PLOT  E 
OR  A1 . IF  P9=6  THEN  NO  PLOT  0 
R  LIST. 

170  '  THE  OUTPUT  ISP, A, W, II,  WH 
ERE  A  CONTAINS  THE  AREA  OF  E 
,  HI  CONTAINS  THE  WINOING  NU 
MBER  U 

175  !  OF  E,  ANO  II  IS  AN  ERROR  I 
NDICATOR .  IF  11=0  OR  2  THEN 
THE  OUTPUT  IS  0.K.  IF  I1=10R 
- 1 , THEN 


180  ‘  ANGULAR  REGION  A1  IS  NOT 
WELL-DEFINED.  IF  I 1=1, THE  VE 
RTEX  OF  R1  AND  ONE  OF  THE  OT 
HER  TWO 

185  *  POINTS  ARE  TOO  CLOSE  TO  EA 
CH  OTHER. IF  I 1 =— 1  THEN  MEASU 
RE  OF  A1  IS  CLOSE  TO  0  OR  2P 
I  . 

190  »  IF  11=2  THEN  2  SIDES  OF  E 
OVERLAP. IF  11=3  THEN  C*C>=1 . 

195  !  IF  11=1  OR  3  THEN  P  IS  SET 
TO  -5.  IF  OUTPUT  IS  O.K.,THE 
N  P  IS  GIVEN  TO  9-DECIMAL  AC 
CURACY. 

208  »  SOURCES :  NSWC^DL  REPORT  #3 
886,  SEPT. 1973.  NSWC'DL  REPO 
RT#80-166,  JUNE  1980.  SIAM  J 
N  SCI. 

205  !  STAT.  COMPUT. , JUNE  1980,  P 
P.179,186.  SIAM  JN.  SCI. STAT 
.  COMPUT, DEC.  1982,  PP.?. 

218  !  X<J)  AND  Y< J)  ARE  DIMENSIO 
NED  AT  90  IF  MORE  POINTS  AR 
E  NEEDED  TO  SPECIFY  E  THEN  M 
AKE 

215  »  CHANGES  AT  LINE  225. 

220  OPTION  BASE  1 

225  DIM  X<90> , YC98) , U1 < 15) 

230  SHORT  19,  J 

235  A 1=2 .71E-19  0  A2= . 90000134  0 
A4=  0007311  0  R2=19. 201924 

240  Rl=. 56418958355  0  R3=4 . 9E-27 
0  T7=PI  0  R4=. 159154943892 
0  T9=l  14472988585  0  R6=.282 
09479177 

245  T8= . 000090000007  8  A7=3 . 2625 
E-l 1  0  Z9=l  12837916709 

250  U1<1>=. 008000408336  0  Ul<2>= 
-.800009718649  0  U1(3>=1  057 
875745E-4  0  U1 <4)=-7 . 0426024 
33E-4 

255  Ul<5>=  003249445432  0  Ul<6>= 
-  811232353215  0  Ul<7>=.0309 
19929552  0  U1 <8>=- . 071490983 
78 

260  Ul< 9 >*.145060043403  0  Ul<10? 
=-.265638206366  0  U1<11>= .44 
2851899329  0  U1 < 12>=- . 666626 
670511 

265  Ul<13)=  886223733187  §  Ul<14 
>=-  999999899776  0  Ul<15)=.8 
36226924931 

278  READ  P8 , P9 , Ml , M2, C, SI , S2, N 

275  C7=1-C*C  0  IF  C7> 0  THEN  285 
ELSE  11=3  0  P=-5 
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280  PRINT  "Il="; II ; "P=";P  8  RETU 
RN 

235  PRINT  " INPUT : “ ; P3; P9; HI >  N2; C 
;S1;S2;N  8  PRINT  8  C7=SQR<C7 

) 

29C  IF  Ntl  THEN  K=N  ELSE  K=3  8  U 
1=0 

295  FOR  J=1  TO  K 

306  IF  P8=0  THEN  READ  X(J>,Y(J> 

305  IF  P9>=1  THEN  PRINT 

>  =  -;X<J>;TBB(16>;  "YC;  J;  ">  =  " 

>  Y(J> 

310  NEXT  J 

315  PRINT  8  IF  P9> 1  OR  P9<0  THEN 
779 

320  FOK  J=1  TO  K 

325  Y< J>=<Y<J>-H2)^S2  8  X(J>=((X 
<J>-H1)/'S1-C*Y(  J>>xC7 

•J7C  ucvT  I 

335  X(K*1>=X<1>  8  Y(K+1 )=Y< 1 ) 

340  P=6  8  11=0  8  0=0  8  K 1=0  8  K= 
1 

345  IF  Nil  THEN  360 
358  X1=X<1)  8  X(3>=X(1>+X<1>-X<3 
)  8  Y<3)=Y< 1 )+Y< 1 >— Y<3)  8  U= 
X<2)-X< 1 )  8  V=Y(2>-Y<1> 

355  W=X< 1 )-X<3)  8  Y1=Y( 1 >  8  Z*Y< 
1 >— Y<3  >  8  GOTO  380 
360  Yl=Y(N-*-l>  8  X1=X(N>1> 

365  U=X(2>—X< 1 >  8  V=Y<2>— Y< 1 ) 

370  X1=X<1>  8  Y1=Y< 1 ) 

375  M=X< 1 >— X(H>  8  Z=Y(1>-Y(N> 

380  Dl=M*H-»-Z*Z 
385  IF  D1>R3  THEN  400 
390  IF  N=1  THEN  730 
395  M=N-1  8  IF  N=2  THEN  665  ELSE 
375 

400  D2=U*U+V*V 
405  IF  02>R3  THEN  430 
410  IF  N=1  THEN  730 
415  K=K+1  8  U=X(K+1>-X1  8  V*Y(K+ 
1 >— Y1  8  D2=U*U+V*V 
420  IF  02<=R3  THEN  415 
425  IF  K=N-1  THEN  665 
430  R=X1*(Y(K+1 >-Y(N>>  8  81*SSR< 
2*01 >  8  B2=SQR<2*02) 

435  P2=V*M-U*Z  8  C1»U*W+V*Z  8  05 
=0TN2(P2,C1>  8  K1=K1+R5  8  L= 
0  8  B= . 5*<X<K)*X<K>*Y<K>*Y<K 
)> 

440  IF  B>01  THEN  450  ELSE  C2»0 
445  P1»05*R4-C2  8  GOTO  645 
450  G1  =  <H*X<K>*Z*Y<K>)/'B1  8  G2»< 
U*X < K > ♦V* Y < K ) > ^B2  8  H1»<Z*X< 
K)-H*Y<K>>/'B1  8  H2=<V*X<K)-U 
*Y<K))/B2 


455  IF  0BS<P2>>2*B1*B2*07  THEN  5 
18 

460  IF  C1<0  THEN  480 

465  IF  0BS<05>  <=T8  THEN  475 

470  IF  G1>=0  THEN  475  ELSE  510 

475  P1=0  8  GOTO  645 

480  11=2 

435  IF  P2<0  THEN  500 
490  H=H2  8  GOSUB  735 
495  Pl= . 5*E1  8  GOTO  645 
508  H=H1  8  GOSUB  735 
595  PI =-< . 5*E1 >  8  GOTO  645 
510  IF  B<=02  THEN  605 
515  IF  G1<0  THEN  545 
520  IF  G2>=0  THEN  625 
525  G2=-G2  8  H2=-H2  8  IF  0BS<H2> 
<=04  THEN  535 

530  H=— H2  8  GOSUB  735  8  L=.5*E1 
8  GOTO  610 

535  L=.5+R1*H2  8  GOTO  618 
548  L= . 5— R1*H1  8  GOTO  610 
545  G1=-G1  8  Hl=— HI 
550  IF  G2<0  THEN  565 
555  IF  BBS (HI >  <=04  THEN  540 
560  H=H1  8  GOSUB  735  8  L=.5*E1  8 
G010  610 

565  G2=-G2  8  H2=~H 2 
570  IF  BBS (HI >  <=04  THEN  590 
575  IF  0BS(H2><=04  THEN  585 
580  H=H1  8  GOSUB  735  8  L=.5*E1  8 
H=H2  8  GOSUB  735  8  L=L- . 5*E 
1  8  GOTO  625 

585  H=H1  8  GOSUB  735  8  L=R1*H2-. 

5*( 1-El )  8  GOTO  625 
598  IF  0BS(H2><=04  THEN  600 
595  H=H2  8  GOSUB  735  8  L=.5*(l-E 
1  )-Rl*Hl  8  GOTO  625 
600  L=R1*(H2-H1>  8  GOTO  625 
605  C2=T6*(H2-H1 )-R4*(G2*H2-Gl*H 
1>  8  GOTO  445 

610  P2*-P2  8  IF  P2<=0  THEN  620 
615  L*L-1  8  05=T7+05  8  GOTO  625 
620  05*05— T7 
625  IF  B>*R2  THEN  640 
630  C3=05  8  C4* . 5*05  8  «*15  8  F* 
0  8  06=K2-H1  8  C5=06  8  GOTO 
715 

635  Pl=L-*-EXP<-<B^T9))*(C4-S3>  8 
GOTO  645 
640  P1*L 

645  IF  KiN  THEN  695 
650  IF  Nil  THEN  675  ELSE  H1=0 
655  IF  K1<*0  THEN  P*0BS(P1)  ELSE 
P*0BS(1-0BS(P1)) 
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660  IF  ABS<K1X. 000000009005  THE 
H  1 1  =- 1 

665  IF  Nil  THEN  PRINT  P;H1;M1;I1 
ELSE  PRINT  Pill 

670  BEEP  0  IF  P9=8  OR  P9=l  THEN 
RETURN  ELSE  I9=P  0  J=H1  6  GO 
TO  1055 

675  P=P-P1  0  K1=K1*R4  0  A= . 5*A  0 
IF  K1 <0  THEN  685 
680  Wl=IP<Kl+.l)  0  GOTO  690 
685  H1=IP<K1-. 1) 

690  P=P+H1  0  H1=S1*S2*C7*A  0  GOT 
0  665 

695  W=U  0  Z=V  @  B 1 =B2  0  X1=X<K+1 
)  0  Yl=Y<k+l)  0  Y2=Y<k> 

700  k=k+l  0  U=X<k+l)-Xl  0  V=Y<k+ 
1)-Y1  0  D2=U*U+V*V 
705  IF  02<=R3  THEN  700 
718  B2=SQR<2*02>  0  P=P-P1  0  A=A+ 
Xl*<Y<k+l)-Y2)  0  GOTO  435 
715  S3=U1<H)*A6 

720  H=H-1  0  H2=H2*G2  0  H1=H1*G1 
0  T=H2-H 1  «  F=F+B  0  C6=<F*C3 
+T)/<16-H>  0  S3=S3+U1 (H)4C6 
725  IF  H=1  THEN  635  ELSE  C3=C5  0 
C5=C6  0  GOTO  720 
730  P=— 5  0  11=1  0  GOTO  665 
735  E1=0  0  C4=ABS<H)  0  IF  HBS<H) 
>=4.4  THEN  755  !  E1=ERFC<H)- 
-9-OECIHAL  DIGIT740  FOR  19=1 
TO  15 

740  E1=<<<<<<<U1<1)*C4+U1<2))*C4 
+U1<3))*C4+U1<4))*C4+U1<5))* 
C4  +UI <6)  )*C4+U1  <7)  )  *C4+U1  <  8) 
>*C4 

745  E1=<<<<<<<E1+U1<9))*C4+U1<10 
) >$C4+U1 < 1 1 >  )*C4+U1 < 12) )*C4+ 
U1 < 13) )*C4+U1 < 14) XC4+U1 < 15) 
)*Z9 

750  E1=E1*EXP<-<H*H>) 

755  IF  H<0  THEN  E1=2-E1 
760  RETURN 

765  BEEP  84/ 100  6  IF  P9=0  THEN  R 
ETURN 

770  U=X<1)  0  V=U  0  W=Y< 1 )  0  2=H 

775  FOR  1=1  TO  K 

780  IF  X<I)>U  THEN  U=X<I) 

785  IF  X<IXV  THEN  V=X<I) 

788  IF  Y(I)>H  THEN  W=YCI> 

795  IF  Y< I ) <2  THEN  2=Y<I) 

890  NEXT  I 

805  GCLEflR  0  P1*U-V  0  IF  Pli0  TH 
EN  810  ELSE  Pl*l  0  V=V-P1  0 
U=U+P1 


810  Y1=LGT <P1 )  0  X1  =  INT < Y1 )  0  Y1 
=FP  <  Y1 ) 

815  Y1=10-Y1  0  IF  XI <0  THEN  Yl=l 
0*Y1 

820  X1=10AX1  0  IF  Yl>=2  THEN  830 
825  Tl= . 1 tXl  0  GOTO  855 
830  IF  Yl>=4  THEN  840 
835  Tl= . 24X1  0  GOTO  855 
840  IF  Yl>  =5  THEN  850 
845  Tl=  4*X1  0  GOTO  855 
850  T1=.5*X1 
855  T3=I NT  <  V^T 1 >  *T 1 
860  T4=U  0  IF  FP<T4/’T1  )#0  THEN  T 
4=<INT<T4^T1)+1)*T1 
865  P2=W-Z  0  IF  P2#0  THEN  870  EL 
SE  P2=l  0  2=2-92  0  W=M+P2 
878  Y1=LGT<P2)  0  X1=INT(Y1)  0  Y1 
=FP< Y1 )  0  Y1=10~Y1  0  IF  XI <0 
THEN  Y1=10*Y1 

875  X1=10^X1  0  IF  Yl>=2  THEN  885 
880  T2=  1*X1  0  GOTO  915 
885  IF  Yl>=4  THEN  900 
890  T2=.2*X1 
895  GOTO  915 
900  IF  Yl>=5  THEN  910 
905  T2=.4*X1  0  GOTO  915 
910  T2= . 5*X1 
915  T5=INT(2/T2)*T2 
920  T6=0+W  0  IF  FP<T6^T2)#0  THEN 
T6= < I NT< T6/T2 ) + 1 ) * T2 
925  SCALE  T3- .  35*P1 ,  T4+- .  25*P1 ,  T5 
-  25*P2/T6-r.2SP2 

930  XAXIS  T5/T1 #T3»T4  0  YAXIS  T3 
/T2/T5/T6 

935  XAXIS  T6,T1,T3.T4  @  YAXIS  T4 
/T2/T5/T6 

940  Gl=T3/<5*Tl>  0  IF  FP<G1)#0  T 
HEN  Gl  =  <INT<Gl)-+l)*5*Tl  ELSE 
G 1 =T3 
945  I =-5 

950  1  =  1+5  0  X=G1+I*T1 

955  IF  X> T4  THEN  980 

960  HOVE  X/T5  0  IORAM  0,.045*P2 

965  HOVE  X»T6  0  I DRAW  0,-<045*P 


978  Y1=LEN<VAL*<X>> 

975  HOVE  X- . 03t < Yl-1 )  tPl , T5- . 1 tP 
2  0  LABEL  VAL*<X)  8  GOTO  958 
980  G1=T5/<5*T2)  0  IF  FP< G1)#0  T 
HEN  G1=(INT<G1 )+l >*5*T2  ELSE 
G1=T5 


983  I=-5 


990  1=1+5  0  Y1=G1+I*T2  0  IF  Y1>T 
6  THEN  1015 

995  HOVE  T3/Y1  0  IDRAW  . 045*P1,0 
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1000  HOVE  T4,Y1  8  I  DRAW  -<045*P 

1  ),0 

1005  X=LEN<VAL*<Y1>> 

1010  HOVE  T3-.05*<l+X>*P3,Yl-.5* 
T2  0  LABEL  VAL*<Y1)  0  GOTO 
996 

1915  PENUP  0  IF  H=1  THEN  1070 
1028  FOR  1=1  TO  N  0  PLOT  X<I>,Y< 
I> 

1025  LABEL  VAL*<I> 

1030  NEXT  I  0  PLOT  X<1>,YC1> 

1035  HOVE  U+l .5»T1,Z-.02*P2  0  LA 
BEL  “X"  0  HOVE  T3- . 02TP1 > T6 
♦  085*P2  0  LABEL  "Y" 

1048  HOVE  T3+.89*P1,T6+  1*P2 
1045  IF  N#1  THEN  LABEL  “SPECIFIE 
D  POLYGON"  ELSE  LABEL  "SPEC 
IFIED  ANG.REG. " 


1050  C8=T3  0  C9=T5-.22*P2  0  D3=P 
1  0  GOTO  320 

1055  HOVE  C8-.2*D3,C9  0  LABEL  "P 
=“fcVAL$<I9>  0  HOVE  C8+.45*D 
3,C9  0  LABEL  "A="fcVAL*<J> 

1060  HOVE  C8+D3,C9  0  LABEL  "W="& 
VAL4CW1 >  0  IF  P9>1  THEN  COP 
Y 

1055  RETURN 

1070  PLOT  X<3>,Y<3>  0  LABEL  VAL* 
<3>  0  PLOT  XO><  Y(l)  0  LABE 
L  VAL*<1> 

1875  PLOT  X<2> » Y<2>  0  LABEL  VALS 
<2>  0  GOTO  1035 

1080  END 
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MLEQRE 

Input:  (Initial  Data  Statement),  N,  M,  LO,  AO,  BO,  P8 

N  :  Number  of  listed  successes 
M  :  Number  of  iisted  failures 

LO  :  If  LO  =  1,  initial  data  statement  is  followed  by  data  statements  containing  list  of 
successes,  aj,  followed  by  data  statements  containing  list  of  failures,  bj,  i.e.,  aj, 
a2>  ....  aN,  b1(  b2,  bM. 

If  LO  ^  1,  initial  data  statement  is  followed  by  data  statements  containing  each 
listed  (different)  success  a*,  with  each  aj  preceded  by  q,  the  number  of  times  that 
aj  occurs.  Then  follow  data  statements  containing  each  listed  (different)  bj  preceded 
by  dj,  the  number  of  times  bj  occurs,  i.e.,  q,  aj,  c2,  a2,  ...,  cN,  aN,  d1?  b1(  d2, 
b2,  ...,  dM,  bM. 

The  listed  successes  (failures)  are  stored  by  MLEQRE  in  Array  A(l)  (S(J)),  I(J)  =  1 
to  N(M).  If  LO  =£  1,  then  the  q(dj)  are  stored  in  Array  C(I)  (D(J)). 

If  LO  =  1  then  ones  are  stored  in  arrays  C(I),  D(J)  by  MLEQRE. 

AO,  BO :  Contain  initial  estimates  for  a  =  p/a  and  0  =  l/o  if  supplied  by  user.  If  BO  <  0, 
then  initial  estimates  for  a  and  0  are  supplied  by  MLEQRE. 

P8  :  If  P8  =  0,  no  confidence  plots  are  made.  If  P8  =  1,  then  confidence  plots  at  the  50 
and  95%  level  appear  on  the  CRT.  If  P8  >  2,  then  plots  appear  on  CRT  and  also  on 
the  printer. 

It  is  necessary  and  sufficient  for  a  solution  that  two  conditions  be  satisfied,  namely 

max  (bj)  >  min  {aj) 
j  i 

(average  of  the  aj)  >  (average  of  the  bj) 


If  either  condition  is  violated  an  exit  is  made  and  a  message  is  printed  indicating  which  condition 
is  violated. 

The  program  is  set  to  limit  the  number  of  iterations  to  30.  This  number  of  iterations  has  never 
been  necessary,  however  if  it  should  occur  the  constraint  can  be  overridden  as  described  in  the 
comment  statements  at  the  head  of  MLEQRE. 

Output:  N,  M.  LO,  AO,  BO,  P8 


//(LO  =  1,  (See  Example  2). 


Successes 

A(  1)  *  3 j 
A(2)  «  a2 

A(N)  *  a. 


Failures 


BO) 

B(2) 

B(M) 


b 

b 


1 

2 

M 
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//( LO  =£  1),  (See  Example  1). 

NO.  OF  A  SUCCESS  AND  SUCCESSES  |  NO.  OF  A  FAILURE  AND  FAILURES 

dy  B(l)  =  by 

d2  B(2)  =  b2 

dM  B(M)  =  bM 

Maximum  likelihood  estimates  MU  and  SIGMA  (SIG). 

Covariance  elements 

auu,  a“»,  a“ 

Final  values  AO  =  n/a,  BO  =  1  la. 

Last  Newton-R?phson  corrections -(contained  in  D1  and  D2). 

Initial  values  for  AO,  BO-(contained  in  A1  and  Bl). 

No.  of  N-R  iterations. 

This  completes  output  if  P8  =  0.  If  P8  =  1  or  >2  plot  follows  on  CRT,  and  in 
the  latter  case  it  also  appears  on  the  printer. 

Accuracy:  Approximately  6-decimal  digits  in  n  and  a  assuming  the  elements  of  (aj,  (bj)  are 
exact. 


Cj  A(I)  =  aj 

c2  A(2)  =  a2 

CN  A(N)  *  a^ 
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N=  4  ,  M=  4  , L0=  0  , A0  =  0  , B0=  .04  , P8=  2 

NC.  OF  ft  SUCCESS  AND  SUCCESSES  NO.  OF  A  FAILURE  AND  FAILURES 

2  A<  1  >*  41.67  6  B  <  1  >=  27.1 

11  A<  2  )  =  58.33  9  B<  2  >  =  41.67 

9  A <  3  ) *  81.67,  10  B (  3  >  =  58.33 

6  A<  4  >*  114.33  1  B<  4  >=  81.67 

MU*  5 . 80059E+0 1  SIG*  1.68263E+01 

COVARIANCE  MATRIX  ELEMENTS 
1 . 28405E+01  ,  1 . 85908E+00  ,  1.99470E+01 

FINAL  A0=  3. 44734E+00  FINAL  B0*  5.94308E-02 
LAST  DELTA  A0,B0  *  3.79341E-08  ,  6.78833E-10 
INITIAL  A0=  0 . 00000E+00  INITIAL  B0*  4.00000E-02 
NO.  OF  ITERATIONS*  7 


MU-5.80059E+01  SIGMR- 1 . 68263E+0 1 

c  50Z  &  95*  CONFIDENCE  ELLIPSES 


N=  10  ,M=  10  ,L0=  1 
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,00=  0  , B3=  0  , F3=  2 


SUCCESSES 
A<  1  >-  2854 
A<  2  >*  2836 
A<  3  >  =  276? 
A<  4  >*  2814 
A<  5  >*  2801 
A<  6  >  =  2792 
A<  7  >-  2820 
A<  8  >*  2741 
A<  9  >*  2767 
A<  10  >*  2761 


FAILURES 
B (  1  >  =  2652 
B <  2  )=  2741 
B<  3  >-  2846 
B<  4  )=  2713 
B <  5  >  =  2806 
B<  6  >  =  2770 
B<  7  >  =  2776 
BC  8  )=  2763 
B<  9  >*  2706 
S (  10  >=  2735 


MU*  2. 77534E+03  SIG*  7.12326E+3* 

COVARIANCE  MATRIX  ELEMENTS 
4 . 63975E+02  , -7 . 89329E+00  ,  1 . 33797E+03 


FINAL  A0=  3 . 896 1 7E+0 1  FINAL  B0=  1.40385E-02 

LAST  DELTA  A0,B0  =  1 . 97272E-05  ,  7.07111E-09 
INITIRL  A0*  5 . 62895E+0 1  INITIAL  B0=  2.02988E-02 

NO.  OF  ITERATIONS-  4 


0 


MU-2.77534E+03  SIGMfl-7. 1232GE+01 

50%  &  95*  CONFIDENCE  ELLIPSES 
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N=  4  ,  M=  6  •  L0=  1  ,  A6=  9  66  = 

>  P8=  2 


SUCCESSES 
A<  1  >=  2481 
fiK  2  >=  2566 
FK  3  >=  2533 
A<.  4  >=  2620 

FAILURES 
B<  1  >=  2443 
B<  2  >=  2486 
B<  3  >=  2463 
B(  4  >=  2480 
B(  5  >=  2505 
B(  6  '=  2500 


MU=2  50393+3  SIG=2 . 62261+i 

COVARIANCE  MATRIX  ELEMENTS 
202  594767688  107  937942839 

338  540527278 


INITIAL  A0=  54  5879207362 
INITIAL  B8=  2. 1 7720294092E-2 
FINAL  80,60=  95  4744043317 
3.81 298920267E-2 
LAST  DELTA  A0,B0= 

1  684434 1 805E-6 
6  74520352426E-10 
N0.  OF  ITERATIONS=  5 
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100  \  THIS  SUBROUTINE,  MLEQRE ,  GIVES  THE  MAXIMUM  LIKELIHOOD  EST¬ 

IMATES, MU  AND  SIGMA,  BASED  ON  OUANTAL  RESPONSE  EXPERIMENTS, 

105  !  POR  THE  MEAN  AND  STANDARD  DEV  I  AT  I  ON , RESPECT  I VELY,  OF  A  NORMAL 

DISTRIBUTION.  OUTPUT  ALSO  INCLUDES  THE  COVARIANCE  ELEMENTS 
110  1  AND  , AT  THE  USER'S  OPTION,  A  PLOT  OF  THE  CONFIDENCE  ELLIPSES 

AT  THE  50  AND  95V.  LEVELS.  THE  QUANTAL  RESPONSES  ARE  CLASSIFIED 
115  1  AS  "SUCCESSES"  OR  "FAILURES".  THE  INPUT  QUANTITIES  TO  MLEQRE 

RESULTING  IN  SUCCESSES( FAILURES)  ARE  ALSO  SIMPLY  REFERRED  TO 
120  1  AS  SUCCESSES  (FAILURES). 

125  !  THE  INPUT  IS  STORED  IN  A  DATA  STATEMENT  BEGINNING  WITH: 

N,M,L0,A0,B0,P8  WHERE  N(M )  DENOTES  THE  HUMBER  OF  "LISTED" 

130  !■  SUCCESSES ( FA  I  LURES ) .  IF  L0*1,  THEN  THE  DATA  STATEMENT  IS 

CONTINUED  WITH  THE  LIST  OF  SUCCESSES  FOLLOWED  BY  THE  LIST 
135  !  OF  "AILURES.  IF  L0#1,  THEN  EACH  DIFFERENT  SUCCESS  IS  PRECEDED 

BY  HE  NO.  OF  TIMES  IT  OCCURS.  THE  SAME  IS  THEN  DONE  WITH 
140  !  FAILURES.  THE  LISTED  SUCCESSES(FAILURES)  ARE  STORED  BY  MLEQRE 

IN  A<IXB<J)>,  I  <  J )  =  1  TO  NCM).  IF  L0#1  THEN  THE  NO.  OF  A(I) 

145  !  ( B  <  J )  >  AT  A  FIXED  I(J)  FOR  EACH  I(J)  IS  STORED  BY  MLEQRE  IN 

C(I)(D(J)>.  IF  L6*  1 ,  THEN  C<IXD<J>>  CONTAINS  N(M)  ONES. 

150  !  A0  AND  B0  CONTAIN  INITIAL  ESTIMATES  FOR  MU^SIGMA 

AND  1/SIGMR.  IF  B0O0,  THEN  A0  AND  B0  ARE  SUPPLIED  BY  MLEQRE. 
155  (MAXIMUM  NO.  OF  DIFFERENT  ACI)  IS  100.  SIMILARLY  FOR  THE  BCJ). 

TO  INCREASE,  CHANGE  DIMENSION  STATEMENT  AT  LINE  345. 

160  !  °8  IS  AN  OUTPUT  SPECIFIER.  IF  P8=0,  THEN  NO  PLOT  IS  MADE.  IF 

P8*l,  THEN  PLOT  APPEARS  ON  CRT.  IF  P8>=2,  THEN. PLOT  APPEARS  ON 
165  !  CRT  AND  PRINTER. 

170  !  IT  IS  NECESSARY  AND  SUFFICIENT  FOR  A  SOLUTION  TO  EXIST  THAT 

MAX.  LISTED  FAILURE>MIN.  LISTED  SUCCESS  AND,  THAT  THE  AVERAGE 
175  !  OF  THE  SUCCESSES >THE  AVERAGE  OF  THE  FAILURES.  IF  EITHER  COND¬ 

ITION  IS  VIOLATED  AN  ERROR  ME.SSAGE  IS  PRINtED  FROM  LINE  710 
180  !  OR  735.  A  MAX.  OF  30  ITERATIONS  IS  ALLOWED.  MESSAGE  PRINTED 

FROM  LINE  800.  IF  MORE  DESIRED  CHANGE'  LINE  790  AND  CONTINUE 
185  !  AT  LINE  770. 

190  !  SOURCES:  NWL  REPORT  2846,  NOV.  1972.  SIAM  JN. 

OF  APPL.  MATH.,  MAR.  1972,  PP.  447-454. 

195  OPTION  BASE  1 

200  DIM  P9<4),Q9<3),R9<6),S9<5),V9<4).W9(3),D*C303,E$C293 

205  P9<i>—  3.56098437018E-2 

210  P9<2>»6. 99638348862 

215  P9<3)«21. 9792616183 

220  P9<4)*242. 667955231 

225  Q9< 1 )-15. 0827976304 

230  Q9<2)«91. 1649054045 

235  Q9<3)»215. 05887587 

240  R9<1).— 6.08581519597E-6 

245  R9<2)». 564371606864 

250  R9<3)*4. 26772010709 

255  R9<4)*14. 5718985969 

260  R9<5)«26. 0947469561 

265  R9<6>»22. 8989928517 

270  S9<l)-7. 56884822936 

275  S9<2)-26. 2887957588 

280  S9<3) 11 50. 2732028633 
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285  S9<4>*51. 9335706876 

290  S9<5>»22. 8989857499 

295  V9( 1 )*-3.24319519278E-2 

300  V9<2)*-. 24391 1029489 

305  V9<3)*-. 119903955268 

310  V9<4>*-. 012130827639 

31b  M9<1)*1. 43771227937 

320  W9<2>*. 489552441961 

325  W9<3>=4. 30026643453E-2 

330  Cl=. 707106781187  ! 1/SQR<2> 

335  C2*. 797884560803  !SQR<2/'P!) 

340  C3=.  564189583548  ISORCl/PI) 

345  DIM  A< 1 00) , B< 1 00 > , C< 100) , DC  1 00) 

350  D*=“NO.  OF  A  SUCCESS  AND  SUCCESSES" 

355  ES=“NO.  OF  fl  FAILURE  AND  FAILURES" 

360  STANDARD 

365  READ  N,  (1,  L0,  A0,  B0,  P8 


370 

8 

375 

PRINT  "N*“;n; “L0*“;L0; ", "A0*“;A0; 

"B0 

PRINT 

380 

IF  L0= 1  THEN  485 

- 

385 

PRINT  D$; TABC38) ; E* 

390 

FOR  1*1  TO  N 

395 

READ  C( I ) , A< I ) 

400 

NEXT  I 

405 

FOR  J*1  TO  M 

410 

READ  DCJ),BCJ) 

415 

NEXT  J 

420 

FOR  1*1  TO  MINCN, M) 

425 

PRINT  C< I ) ; TABC9) ; "A<"; I ; "  >*" ;  AC  I > ; TA3< 38) ; DC  I  > ; 

TAB  <  45  > 

430 

NEXT  I 

435 

IF  M>N  THEN  465 

440 

IF  N*N  THEN  590 

445 

FOR  I *M+ 1  TO  N 

450 

print  C<I);tab<9>;  "A<";  I;  *>*••;  aciJ 

455 

NEXT  I 

460 

GOTO  590 

465 

FOR  I-N+l  TO  M 

470 

PRINT  TAB<38);D<I);TAB<45); mb<h; i; ")»“;B<I> 

475 

NEXT  I 

480 

GOTO  590 

485 

PRINT  D*C22,303;TAB<30);E*t22,293 

490 

FOR  I»,l  TO  N 

495 

READ  A< I ) 

500 

C<  I  >* 1 

505 

NEXT  I 

510 

FOR  J-l  TO  M 

515 

READ  B<  J) 

520 

D<J)-1 

• 

525 

NEXT  J 

530 

FOR  1*1  TO  MINCM.N) 

535 

PRINT  MA<"; I; ")-";A<l);TAB<30>: "B<"; i; ">*":B<I> 

540 

NEXT  I 
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345  IF  M  >N  THEN  575 

550  IF  N*N  THEN  590 

555  FOR  I*MM  TO  N 

560  print  “A<  I ;  •• )-“ ; ft <  I ) 

565  NEXT  I 

570  GOTO  590 

375  FOR  I*N+1  TO  M 

580  PRINT  TRB<30>;  “B<“;  i;  = 

585  NEXT  I 
590  !  PAUSE 

595  Z*1E99 

600  D*-Z 
6e5  El*. 008 
610  E2*. 01 

615  L 1*L2*L3*0 
620  FOR  1*1  TON 
625  L1*L1+C< I ) 

630  IF  flCIXZ  THEN  Z*A<I) 

635  P1*C< I >*R< I > 

640  L2*L2+P 1 

645  P2-Pl*fl<I) 

630  L3-L3+P2 

655  NEXT  I 

660  L4=L5*L6*0 

665  FOR  J*1  TO  M 
670  L4*L4+D< J> 

675  IF  B<J>>D  THEN  D*B<J> 

680  P3*D< J)»B< 

685  L5-L3+P3 

690  P4=P3*B<J> 

695  L6*L6+P4 

700  NEXT  J 

703  IF  2<D  THEN  720  , 

710  PRINT  "R-NJN>*B-MRXM 

713  RETURN 

720  L7*L2/L1 

725  L8*L5^L4 

730  IF  L7>L8  THEN  745 

735  PRINT  "RVE.  OF  fl's  <*  RVE.  OF  B's“ 

740  RETURN 

745  IF  B0>0  THEN  770 

750  L1*L1+L4 

753  V*<L2+L3>'L1 

760  B0*SQR<1/<<L3+L6)/'L1-V*V)) 

765  R0*.3*<L?+L8>*B0 

770  fil*A0 

775  B1*B0 

780  Z-L6*L7*L8«I.9*0 

785  L1*L2*L3*L4*L3*0 

790  IF  L9< 30  THEN  810 
795  U*1E99 

800  PRINT  “30  ITERATIONS* 

805  RETURN 
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810 

L9-L9+1 

813 

FOR  I«1  TO 

N 

820 

S*B0*fl<I 

>-fl0 

823 

T  *-S 

838 

GOSUB  1215 

835 

P*C(I)iV 

840 

Ll-Ll-P 

843 

Pl*fl<l>*P 

850 

L2«L2+P1 

835 

P2*C< I >*Y1 

860 

L3-L3-P2 

865 

P3=fl< I >*P2 

870 

L4-L4+P3 

873 

P4*fl< I >*P3 

880 

L5-L5-P4 

885 

IF  Z«0  THEN  920 

890 

P5*C< I )*Y2 

895 

L6-L6+P5 

900 

P6*S*P5 

905 

L7-L7+P6 

910 

P7-S*P6 

915 

L8-L8+P7 

920 

NEXT  I 

923 

FOR  J»1  TO 

M 

930 

T-B0*B<J 

>-R0 

935 

GOSUB  1215 

940 

P*D< J)*Y 

945 

Ll-Ll+P 

930 

P1«B<J>*P 

953 

L2*L2-P1 

960 

P2*D< J)*Y1 

963 

L3-L3-P2 

970 

P3*B< J)*P2 

973 

L4-L4+P3 

980 

P4-BCJ)* 

>3 

983 

L5-L5-P4 

990 

IF  Z-8  T 

HEN  1025 

993 

P5«D< J>* 

Y2 

1000 

L6-L6+P5 

1005 

P6-T#P3 

1010 

L7-L7+P6 

1015 

P7-T*P6 

1020 

L8-L8+P7 

1025 

NEXT  J 

1030 

D»L3*L5-L4 

*L4 

1035 

D1-<L2*L4- 

LHM-3VD 

1040 

D2“<L1»L4- 

L2*L3>'D 

1045 

R0-A0+D1 

1050 

B0-B0+D2 

1033  IF  <fiBS<Dl > >E1*<PBS<R0>*. 00000001))  OR  <RBS<D2> >E2*B0>  THEN  783 

1060  Z-Z+l 

1063  IF  Z<2  THEpI  783 

1070  S«1'B0  IS^SIGMH 
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1073  U»A0*S  ! U-MU 

1080  !  L6-Auu,  L7*Aus,  L8-A*s 

1085  2-B0*B8 

1090  L6»L6*Z 

1093  L7-L7*2 

1100  L8-L8*Z 

1105  L0-L6*L8-L7*L7 

1110  L1M/-L0 

1113  L2-L8*L1 

1120  L3— L7*L1 

1125  L4-L6*L 1 

1130  !  L2,L3,L4  ARE  THE  COVARIANCE  MATRIX  ELEMENTS. 

1133  PRINT 
1140  FLOAT  ,3 

1145  print  “mu*";u;tab<19>; "Sic*";s 
1150  PRINT 

1155  PRINT  "COVARIANCE  MATRIX  ELEMENTS  “ 

1160  print  L2; ", ";ls; ", ";L4 
1163  PRINT 

1178  PRINT  "FINAL  A0-" ; A0; TAB<26> J “FINAL  B0-";B0 
1175  PRINT  “LAST  DELTA  A0,B0  »“;D1;",";D2 
1180  PRINT  “INITIAL  A0-" ; A1 ; TAB<27) ;" INITIAL  B0-";B1 
1185  STANDARD 

1190  PRINT  "NO.  OF  ITERATIONS-" J L9 
1195  PRINT 

1200  IF  P8 >0  THEN  1440  !  FOR  CONFIDENCE  ELLIPSE  PLOTS. 

1205  RETURN 

1210  !  CODY.  INPUT1T.  OUTPUT!  Y-SQR<2'PI  >*EXP<-Kl»Ki  >/ERFC<Kl  )  ,  Y1*Y*<Y-T>, 

Y2-  SQR<2/PI>*EXP<-Kl*Kl>*Y/<2-ERFC<Kl>>,  12-DECIMAL  ACCURACY. 

1213  K1-T*C1 

1220  Y-Y1-Y2-0 

1223  IF  K 1 < --5. 5  THEN  RETURN 

1230  IF  ABS ( K 1 >  > . 5  THEN  1273 

1235  K4-K 1 *K 1 

1240  K3=<<<P9d)*K4  +  P9<2))*K4-*-P9<3))*K4+P9<4>>/<<<K4+Q9<l>)*K4+Q9<2))»K4+Q9<3)> 

1245  Y- 1 -K 1 *K3 
1250  W-C2*EXP  < -K4  ) 

1253  Y-U/Y 
1260  Y 1 -Y*  <  Y-T  > 

1265  IF  200  THEN  Y2-W*Y/ <  1 +K 1  *K3 > 

1270  RETURN 
1275  K4-A3SCK1) 

1280  IF  K4 >4  THEN  1350 

1283  K3-<  <<<  <R9<  1  >*K4  +  R9<2:>  >*K4+R9<3>  >*K4  +  R9<4>  >*K4+R9<5>)*K4+R9<6>  >✓<<  <  <  CK4  +  S9 

<1))*K4+S9<2>>*K4+S9<3>)*K4+S9<4>)*K4+S9<5>> 

1290  IF  K 1  <  0  THEN  1320  v 

1295  Y-C2/K3 

1300  IF  2-0  THEN  1340 

1303  W-EXP  < -K 1  *K  1  > 

1310  Y2»C2*W*Y/<2-W*K3> 

1315  GOTO  1340 
1320  W-EXP  < -K 1 *K  1  ) 

1323  Y-C2*W^(2-W#K3> 

1330  IF  2-0  THEN  1340 
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1333  Y2«C2*Y/K3 

1340  Y1»Y*<Y-T> 

1343  RETURN 
1350  K6*1^<K1*K1 ) 

1333  K3»<  < (V9< 1 >*K6  +  V9<2>  >*K6+V9<3>  >*K6+V9<4>  >/<  <  <K6+W9< 1 >  >*K6+W9<2>  >*K6+U9<3>  > 

1360  U1*C3+K3*K6 

1363  IF  K1<0  THEN  1400 

1370  Y«C2*K1/W1 

1375  Yl  — <2#C3*K3/'<W1*W1>> 

1380  IF  <Z»0>  OR  <K1 >5. 5)  THEN  RETURN 
1385  W»EXP<-K1#K1 ) 

1390  Y2*C2*U*VV<2-W*W1/K1 > 

1395  RETURN 

1400  U*EXP< -K 1  *  K 1  > 

1405  Y»C2*W/<2“W*W1/'K4J 

1410  Y1«Y*<Y-T> 

1415  IF  Z-0  THEN  RETURN 
1420  Y2«C2*K4#Y/W1 

1425  RETURN 
1430  Y*Y 1 *Y2*0 
1435  RETURN 

1440  EXIT  GRAPHICS  !  PLOTTING  OF  CONFIDENCE  ELLIPSES  FOLLOWS. 

1443  PLOTTER  IS  “GRAPHICS" 

1430  GRAPHICS 

1453  L9-5.99  !L9»  95*  CHI-SQUARED  VALUE. 

1460  S1-2*SQR<L8#L9^L0>  !  RANGE  OF  X-VALUES. 

1463  S2»2*SQR<L6*L9^L0>  !  RANGE  OF  Y-VALUES. 

1470  Y-LGT<S1> 

1473  T-INT<Y) 

1480  Y"FRACT  < Y> 

1485  Y»10‘NY 

1490  IF  T<0  THEN  Y»10#Y 
1493  T-10"^ 

1500  IF  Y>-2  THEN  1313 

1303  K1-.1*T 

1310  GOTO  1330 

1315  IF  Y>4  THEN  1330 

1320  K1».2*T 

1323  GOTO  1530 

1330  IF  Y>»3  THEN  1343 

1333  K1».4#T 

1340  GOTO  1330 

1343  K 1 ■ . 3*T 

1330  K3-INT<<U-.5*S1)/K1>*K1 

1535  K4«U+.3#S1 

1360  IF  FRACT <K4/K1><>0  THEN  K4-< INT<K4/K1 >♦! >*Ki 
1363  Y-LGT<S2> 

1370  T«INT<Y> 

1373  Y»FRACT<Y> 

1380  YM0*Y 

1383  IF  T<0  THEN  Y-10#Y 

1390  T-IO-'T 

1593  IF  Y>«2  THEN  1610 
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1600  K2«.1*T 

1605  GOTO  1645 

1610  IF  Y>»4  THEN  1625 

1615  K2».2*T 

1620  GOTO  1645 

1625  IF  Y>«5  THEN  1640 

1630  K2«.4*T 

1635  GOTO  1645 

1640  K2«.5*T 

1645  K5»INT<<S-.5*S2>/k2>«K2 

1650  K6*S+.5#S2 

1655  IF  FRACT<K6/K2><>0  THEN  K6-<  IN7<K6/K2>+1  )*K2 

1660  SCALE  K3-.35*Sl,K4+,25*Sl,K5-.25*S2,K6+.2*32 

1665  CLIP  K3,K4,K5,K6 

1670  LAXES  K1,K2,K3,K5,5,5,4 

1675  LAXES  K1 , K2, K4, K6, 5, 5, 2 

1680  !  ELLIPSE  PLOT 

1685  Ll-Cl 

1690  L2-SGN<L7)*C1 

1695  PI *L6-L8 

1700  IF  ABS<P1K<ABS<L6)+ABS<L8))*1E-10  THEN  1725 
1705  P3»ATN<2*L7/P1) 

1710  P3».5*P3 

1715  L1«C0S<P3> 

1720  L2*SIN<P3) 

1725  1*0 

1730  K»40 

1735  P5*P I +P I 

1740  P7*P5/'K  !  K*  N0.  OF  POINTS  FOR  ELLIPSE. 

1745  X*L 1 *L 1 
1750  Y*L2*L2 

1755  T»SQR<L9/'<L6*X+L8#Y  +  2*L7*L1*L2>  > 

1760  Y1-SGS<L9-'<L6*Y  +  L8*X-2*L7»L1*L2>  > 

1765  X*COS  <  P7  > 

1770  Y»SIN<P7) 

1775  P3*T#L2 

1780  P4»Y1»L1 

1785  Pl-TirLl 
1790  P2*Y 1 *L2 
1795  Z»-P7 

1800  W* 1 
1805  Y2-0 

1810  MOVE  U+P1.S+P3 
1815  Z-Z+P7 
1820  IF  Z>P5  THEN  1850 
1825  P6»W#X-Y2*Y 

1830  Y2»W*Y+Y2*X 

1835  W*P6 

1840  DRAW  U*P1*W-P2*Y2,S+P3#U+P4*Y2 

1845  GOTO  1815 

1850  IF  IO0  THEN  1900 

1855  I«1 

1860  MOVE  U+Pl , S+P3 
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365  DRAW  U-Pl , S-P3 
370  nuVE  U+P2.S-P4 
375  DRAW  U-P2.S+P4 

380  I8=SQR< 1 . 39^L9)  11.39  IS  50*  CHI-SQUARED  VALUE. 

385  T*T*I8 
390  Y 1 *Y1 *  1 8 
GOTO  1775 

300  MOVE  K3+. 15*S1,K6+.025*S2 

305  LABEL  “50*  «,  95*  CONFIDENCE  ELLIPSES" 

310  FLOAT  5 

315  MOVE  K3-.005*S1,K6+.08#S2 
320  LABEL  “MU»“*.VAL*<U> 

325  MOVE  K4-.4*S1,K6+.08*S2 
330  LABEL  “SIGMA«“8,VAL*<S> 


335 

,  MOVE  K4+ . 02*S 1 , K5 

340 

LABEL  “U“ 

345 

MOVE  K3-.005*S1,K6* 

.  01*S2 

350 

LABEL  "$“ 

355 

IF  P8>-2  THEN  DUMP 

GRAPHICS 

360 

RETURN 
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100  '  THIS  SUBROUTINE,  MLEQRE,  G 
IVES  THE  MAXIMUM  LIKELIHOOD 
ESTIMATES,  MU  AND  SIGMA,  BAS 
ED  ON 

105  »  QUANTAL  RESPONSE  EXPERIMEN 
TS,  FOR  THE  MEAN  AND  STANDAR 
D  DEVIATION, RESPECTIVELY,  OF 
A  NORMAL 

110  !  DISTRIBUTION.  OUTPUT  ALSO 
INCLUDES  THE  COVARIANCE  ELEM 
ENTS  AND  , AT  THE  USER'S  OPTI 
ON,  A  PLOT 

115  »  OF  tHE  CONFIDENCE  ELLIPSES 
AT  THE  50  AND  95*  LEVELS.  T 
HE  QUANTAL  RESPONSES  ARE  CLH 
SSIFIED 

120  »  AS  -SUCCESSES"  OR  "FAILURE 
S“.  THE  INPUT  QUANTITIES  TO 
MLEQRE  RESULTING  FROM  SUCCES 
SES< 

125  •  FAILURES)  ARE  ALSO  SIMPLY 
REFERRED  TO  AS  SUCCESSESCFAI 
LURES).  THE  INPUT  IS  STORED 
IN  A 

130  '  OATA  STATEMENT (S)  BEGINNIN 
G  WITH  N,M,L0, A0,B0,P8,^HERE 
N(M)  DENOTES 

135  !  THE  NUMBER  OF  "LISTED"  SUC 
CESSES(FAILURES) .  IF  L0=i,  T 
HEN  THE  DATA  STATEMENT  IS  CO 
NTIHUED 

140  *  WITH  THE  LIST  OF  SUCCESSES 
FOLLOWED  BY  THE  LIST  OF  FAX 
LURES.  IF  L0#1 ,  THEN  EACH  DI 
FFERENT 

145  •  SUCCESS  IS  PRECEDED  BY  THE 
NO.  OF  TIMES  IT  OCCURS.  THE 
SAME  IS  THEN  DONE  WITH  THE 
FAILURES. 

150  !  THE  LISTED  SUCCESSES < FA I LU 
RES)  ARE  STORED  BY  MLEQRE  IN 
ACIXBCJ)),  I<J)*1  TO  N<M>  . 
IF  L0t 1 

155  ?  THEN  THE  NO.  OF  fl(IXB(J)) 
AT  A  FIXED  I<J>  FOR  EACH  I< 
J>  IS  STORED  BY  MLEQRE  IN  C< 
I ) <D< J> ) 

160  •  IF  L0= 1 ,  THEN  C<IXn<J>>  C 
ONTAINS  H<M >  ONES. 

165  f  A0  ANO  B0  CONTAIN  THE  INIT 
IAL  ESTIMATES  FOR  HU'SIGHA  A 
HD  1/SIGMA,  RESPECTIVELY 

170  !  IF  B0<»0,  THEN  MLEQRE  SUPP 
LIES  THE  INITIAL  ESTIMATES  F 
OR  A0  AND  B0 . 


175  !  THE  A<I),B<J),C<I),D< J>  AR 
E  DIMENSIONED  100  AT  LINE  27 

5. 

180  !  P8  IS  AN  OUTPUT  SPECIFIER. 
IF  P8=0 , THEN  NO  PLOT  IS  MADE 
.  IF  P8=l , THEN  PLOT  APPEARS 
ON  CRT. 

185  !  IF  P8>=2,  THEN  PLOT  APPEAR 
S  ON  CRT  AND  ON  THE  PRINTER. 

190  •  IT  IS  NECESSARY  AND  SUFFIC 
IENT  FOR  A  SOLUTION  TO  EXIST 
THAT  THE  MAX.  FAILURE>MIN 
SUCCESS , 

195  »  AND  THAT  THE  AVERAGE  OF  TH 
E  SUCCESSES) AVERAGE  OF  THE  F 
AILURES .  IF  EITHER  CONDITION 
IS  VIO- 

200  !  LATED  AN  ERROR  MESSAGE  IS 
PRINTED  FROM  LINES  440  OR  45 
0.  A  MAX. OF  30  ITERATIONS  IS 
ALLOWED 

205  •  IF  MORE  DESIRED, CHANGE  30 
IN  LINE  480  AND  CONTINUE  AT 
470. 

213  !  SOURCES  NWL  REPORT  2846, NO 
V . 19720SI AM  JN  OF  APPL .  MAT 
H.,  MAR  1972, PP  447,454. 

215  OPTION  BASE  1 

220  DIM  P9(4> , Q9<3> , R9(6> , S9<5> , 
V9  <  4  > , W9  <  3 ) , DSC30U , E*C293 

225  P9< 1 )=“3 . 560984370 13E-2  0  P9 
<2 >=6.99638348862  0  P9<3>=21 
.9792616183  0  P9< 4 > =242 . 6679 
55231 

230  Q9< 1 >*15 . 0827976304  0  Q9<2>= 
91.1649054045  0  Q9<3)=215.05 
887587 

235  R9<l)*-6  08581519597E-6  0  R9 
<2>= .564371686864  0  R9<3>=4. 
26772010789 

240  R9<4)=J4. 5718985969  0  R9<5>= 
26 . 094746V561  0  R9<6)=22.898 
9928517 

245  S9<1 >*7.56884822936  0  S9<2>= 
26.2887957588 

250  S9<3>=50  2732028638  0  S9<4>= 
51.9335786876  0  S9X5)*22 . 898 
9857499 

255  V9<  1  >**-3 . 24319519278E-2  0  V9 
<2)=- . 24391 102S489  0  V9<3)=- 
.119903955268  0  V9<4)=~.0121 
30827639 

260  W9<1)*1 .43771227937  0  W9<2>* 
.489552441961  0  W9<3>*4.3802 
6643453E-2 
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Cl=. 707106781 187  0  C2= . 79788 
4560803  8  C3=. 564 189583548  8 
•  C2=SQR<2 'PI )0C3=SQR< l^PI ) 

•  INPUT  N,rl,L0>A6,B0 
DIN  A< 100) >C(1 00  >  j  B  < 1 00  )  >  D(  1 
00) 

READ  N,H,L0, A0,B0,P8 
PRINT  "N=";N;*, " ; "H=" ; H; " , " ; 
"L0="; ;L0; ", "A0="; A@; ", ■ 
B0=";B0; ", ";*P8=“,P8  0  PRINT 
D*=“N0.  OF  R  SUCCESS  RND  SUC 
CESSES"  0  £*="N0.  OF  fl  FRILU 
RE  AND  FAILURES" 

IF  L0=  1  THEN  350 
PRINT  D* 

FOP  1=1  TO  N 
READ  C<I),A<I) 

PRINT  C<I)iTAB<9>; "A<"> Ij “>  = 
" ; A< I  ) 

NEXT  I 

PRINT  8  PRINT  E$ 

FOP  J=1  TO  H 
READ  D(J)»B(J) 

PRINT  D<J);TAB<9)i "B<"; J; “)= 

“:B<J)  0  NEXT  J 

GOTO  388 

PRINT  DSC22, 302 

FOR  1=1  TO  N 

READ  A  <  I  )  0  PRINT  "A<";I;")=" 
A  (  I  )  e  C<I>=1  0  NEXT  I 
PRINT  0  PRINT  ESC22, 293 
FOR  J=1  TO  H 

READ  B<J)0  PRINT  "B<";J;">=" 
;B<J>  0  0< J)=l  0  NEXT  J 
Z=1.E99  0  D=-2  0  El=  008  @  E 
2=  01  0  LET  LI , L2,  L3=0 
FOR  1=1  TO  N 

L 1 =L 1 +C <  T  )  0  IF  A< I ) <2  THEN 
2=A< I ) 

P1=C< I >*R< 1 )  0  L2=L2+P1 
P2=P1*A<I>  0  L3=L3+P2 
NEXT  I 

L4=0  0  L5=0  0  L6=0 
FOR  J=1  TO  M 

L4=L4+0<  J)  0  IF  B<J)>D  THEN 
D=8< J> 

P3=0<J)*8<J>  0  L5=L5+P3 
P4=P3*6<J>  0  L6=L6+P4 
NEXT  J 

IF  2>=0  THEN  PRINT  *A-HIN>=B 

-MAX"  ELSE  458 

RETURN 

L7=L2^L  1  0  L8=L5/'L4  0  IF  L7< 
=L8  THEN  PRINT  “RUE  OF  THE  A 
‘ S <=AVE  OF  THE  B’S"  ELSE  460 


455  RETURN 

460  IF  B0>0  THEN  470 
465  L1=L1+L4  0  V= (L2+L5 ) ^L 1  0  B0 
=SQRU/<  CL3+L6)^L1-V*V>)  0  A 
0= . 5*CL?+L8>*B0 
470  A1=A0  0  B 1 =B0 
475  LET  Z, L6 , L7 , L8  , L9=0 
480  LET  L1,L2,L3,L4,L5=0  0  IF  L9 
<30  THEN  L9=L9+1  ELSE  U=1.E9 
9  0  PRINT  "30  ITERATIONS"  0 
RETURN 

485  FOR  1=1  TO  N 
490  S=B0*A<I>-A0  0  T=-S 
495  GOSUB  695 
506  P=C<I)*Y  0  L1=L1-P 
505  P1=A< I  >*P  0  L2=L2+P1 
518  P2=C<I)*Y1  0  L3=L3-P2 
515  P3=A< I > *P2  0  L4=L4+P3 
520  P4=A<I)*P3  0  L5=L5-P4 
525  IF  Z=0  THEN  545 
530  P5=C  < I )*Y2  0  L6=L6+P5 
535  P6=S»P5  0  L7=L7+P6 
540  P?=S*P6  0  L8=L8+P7 
545  NEXT  I 
550  FOR  J=1  TO  H 
555  T=B0*B<J)-A0 
560  GOSUB  695 
565  P=D<J)*Y  0  L1=L1+P 
570  P1=B<J)*P  0  L2=L2— PI 
575  P2=0<J)*Y1  0  L3=L3-P2 
580  P3=B<J)*P2  0  L4=L4+P3 
585  P4=B<J)*P3  0  L5=L5-P4 
590  IF  Z=0  THEN  610 
595  P5=0< J)*Y2  0  L6=L6+P5 
600  P6=T*P5  0  L7=L7+P6 
685  P7=T*P6  0  L8=L8+P7 
618  NEXT  J 

615  0=L3*L5-L4*L4  0  01=<L2*L4-L1 
*L5)/'D  0  D2=<L1*L4-L2*L3)/D 
620  A0=A0+D1  0  B0=B0+O2 
625  IF  ABS<D1 >>E1*<ABS< A0>+ . 8000 
0001)  OR  ABS<D2)>E2*B0  THEN 
480 

638  Z=Z+1  @  IF  Z42  THEN  488 
635  S=1/B0  0  U=A0*S 
64e  *  U=HU@S=S I GNA0L6=A  *  SUB  MUMU 
1 0L7=R ' SUB  HUS ' @L8=A  * SUBSS  * . 
645  Z=B0*B0  0  L6=L6*Z  0  L7=L7*Z 
0  L3=L8*Z  0  L0=L6*L8-L7*L7  0 
L1=1'L8  0  L2=L8*L1  0  L.3=-<L 
7*L1)  0  L4=L6*L1 
653  *  L2, L3,  L4  ARE  THE  ELEMENTS 
OF  THE  COVARIANCE  MATRIX 
655  PRINT  0  Z1=U  0  GOSUB  1160 
666  B*=A*  0  Z1=S  8  GOSUB  1168 
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65  PRINT  "HU=";BS;TAB<17>i "SIG= 
";A$  0  PRINT 

79  PRINT  “COVARIANCE  MATRIX  ELE 
MENTS-  0  PRINT  L2;L3;L4  0  PR 
INT 

75  PRINT  "INITIAL  A0=“;A1 , "IN1T 
IAL  B0=" ; BT > "FINAL  A0,B8=";A 
9,B8,"LHST  DELTA  A0,B9=";D1; 
D2 

580  PRINT  “N0 .  OF  ITERATIONS=" ; L 

9  0  PRINT  0  IF  P8#8  THEN  825 
ELSE  RETURN  !  CONFIDENCE  PL 

OTS. 

>85  REM  CODY.  INPUT  T.  OUTPUT: 

Y<K1 >=SQR<2/PI >tEXP<-KltKl> 
/ERFCCKl),  Yl=Yt<Y-T) 

390  !  OUTPUT  CONTINUED  Y2=SQR<2/ 
PI >*EXP<-K1*K1 >tY/<2-ERFC> 
595  Kl=TtCl  0  Y=0  0  Yl=8  «  Y2=6 
0  IF  K 1 <=— 5 . 5  THEN  RETURN 
700  IF  ABS CK 1 ) > .5  THEN  730 
795  K4=KltKl 

710  K3=<<<P9<l)tK4+P9<2>>tK4+P9< 
3>>tK4-fP9<4>>/<<<K4+Q9<l)  >tK 
4+Q9<2>)tK4-fQ9<3>> 

715  Y=l-KltK3  0  W=C2tEXP<-K4>  0 
Y=U/Y  0  Yl=Yt<Y-T> 

720  IF  Z#0  THEN  Y2=UtY^< 1+K1 tK3> 
725  RETURN 

730  K4=A8S<K1)  0  IF  K4>4  THEN  77 
5 

735  K3=<  < <  <R9< 1 >tK4+R9<2> > tK4+R9 
<3> )*K4+R9<4) )*K4+R9<5) >tK4+ 
R9<6> 

748  K3=K3/<  <  < <  <K4+S9< 1 >  >tK4+S9<2 
>  >tK4+S9<3> )*K4+S9<4) >tK4+S9 
<. 5>  ) 

745  IF  K 1  >  =0  THEN  Y=C2/'K3  ELSE  7 
60 

759  IF  Z#0  THEN  W=EXP<-<K1*K1 > > 
ELSE  778 

755  Y2-C.tUtY/<2-Wtk3>  0  GOTO  77 
0 

768  W=EXP<-<KltKl>>  0  Y=C2*M/<2- 
WtK3  > 

765  IF  Z#0  THEN  Y2=C2tY/K3 
770  Yl=Yt<Y-T>  0  RETURN 
775  K6*l/<KltKl) 

780  K3=<  <  (  V9 ’>1  )tK6+V9<2)  ) tK6+ V9< 
3) >*K6+V9<4>>/<  <<K6+H9< 1 >>tK 
6+W9<2> >*K6+W9<3> ) 

785  Ul=C3+K6tK3  0  IF  KJ <0  THEN  8 

10 

799  Y=C2*Kl/i41  0  Yl=-<2tC3tK3/< W 

ItWli) 


795  IF  ?.#0  THEN  W=EXP<-<KltKl > > 
ELSE  RETURN 

808  Y2=C2*W*Y/<2-W*Ml'Kl>  0  RETU 
RN 

805  Y=0  0  Y 1 =0  0  Y2=0  0  RETURN 
819  W=EXP<-aa*Kl>>  0  Y=C2*W/<2- 
UtWl/K4>  C  Yl=Yt<Y-T> 

815  IF  Z#0  THEN  Y2=YtC2tK4/'Ul 

829  RETURN 

825  ‘  PLOTTING  BEGINS 

830  GCLEAR  §  L9=5 . 99  »  L9=95*:CHI 
-SQUARED. 

835  S1=2*SQR<L3*L9/L0>  0  S2=2tSQ 
R<L6*L9' L0  > 

840  Y=LGT<S1V  0  T=INTCY>  0  Y=FP< 

Y  > 

345  IF  T<(2  Y=10A(Y+1)  ELSE 

Y=10-Y 
850  T=10-T 

855  IF  Y>=2  THEN  865 

869  Kl=.ltT  0  GOTO  399 
865  IF  Y>=4  THEN  875 

870  Kl=.2tT  0  GOTO  899 
875  IF  Y>-5  THEN  885 

889  Kl=.4tT  0  GOTO  899 
885  Kl=.5tT 

890  K3=INT^(U-.5*S1>/K1>*K1 
895  K4=U+.5tSl  0  IF  FP<K4/K1>#0 

THEN  K4=(INT(K4^K1>+1>*K1 
906  Y=LGT(S2)  0  T=INT< Y>  0  Y=FP( 
Y) 

905  IF  T <0  THEN  Y=10^<YH>  ELSE 
Y=10~Y 
910  T=10-T 

915  IF  Y>  =  2  THEN  925 
929  K2=.ltT  0  GOTO  950 
925  IF  Y>=4  THEN  935 
936  K2=.2tT  0  GOTO  958 
935  IF  Y>=5  THEN  945 
949  K2=.4tT  0  GOTO  958 
945  K2=.5*T 

359  K5  =  INT<  <S- .  5»S2)/'K2)*K2 
955  K6=S+.5tS2  0  IF  FP<K6/K2>#0 
THEN  K6=<INT<K6/K2>+1)*K2 
968  SCALE  K3- . 35tSl , K4+ . 25tSl , K5 
-.25tS2.K6-f.2tS2 

965  XAXIS  K5.K1.K3.K4  0  YAXIS  K3 
.K2.K5.K6 

970  XAXIS  K6.K1.K3.K4  0  YAXIS  K4 
.K2.K5.K6 

975  Wl=K3/(5tKl>  0  IF  FP< W1 >#8  T 
KEN  Wl  =  < INT <W1 )  +  l ) tStKl  ELSE 
N 1  =K3 
980  J=-5 
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J= J+5  0  X=W1+J*K1  0  IF  X> K4 
THEN  1010 

MOVE  X,K5  0  I DRAW  0,.045*S2 
MOVE  X,K6  0  IDRftW  0,-<045*S 
2> 

1^Y=LEN*;  VAL$<X> ) 

.  MOVE  X-.03*<Y-1>*S1,K5-  ITS 
2  0  LftBEL  VALS<X>  0  GOTO  98 
5 

i  W1=K5/<5TK2>  0  IF  FPCW1>#9 
THEN  Wl  =  < I NT €W1 )+l )T5TK2  EL 
SE  W1=K5 
J=-5 

i  J=J+5  0  Y=W1 +JTK2  0  IF  Y>K6 
THEN  1050 

MOVE  K'3 /  Y  0  IDRftW  045TS1.0 
MOVE  K4,Y  0  IDRAL  -<845*S1 

>,0 

X=LEN<  VftL$  <  Y  >  > 

MOVE  K3-  05»U+X>TS1,Y-.5*K 
2  0  LftBEL  VflLS1' Y)  0  GOTO  19 
20 

•  ELLIPSE  PLOTS 
L1=C1  0  L2=SGN<L7)*C1  0  P 1= 
L6-L8  0  IF  ftBS <Pl><<ftBS<L6> 
+flBS<L8>>*. 0000090001  THEN 
I960 

P3= . 5TflTN<2*L7/Pl )  0  Ll=COS 
<P3>  0  L2=S I N < P3 > 

1=0  0  K=30  0  P5=PI^PI  0  P7= 
P5/K  !  K=  N0 .  OF  POINTS  FOR 
PLOTS  OF  ELLIPSIS. 

X=L 1 TL 1  0  Y=L2*L2  0  T=SQR<L 
9'( L6*X+L8*Y+2*L7TL1 TL2) >  0 
Y1=SQR<L9^<L6*Y+L8TX-2TL7T 
L 1 TL2 ) > 

X=COS<P7>  0  Y=S I N< P7 ) 

P3=T*L2  0  P4=Y1*L1  0  P1=T*L 
l  0  P2=Y1*L2 
2=-P7  0  W=1  0  Y2=0 
MOVE  U+P1,S+P3 
2=Z+P?  0  IF  2>P5  THEN  1105 
P6=M*X-Y2*Y  0  Y2=W*Y+Y2*X  0  ~ 
M=P6  0  DRHW  U-*-Pl*U-P2*Y2,S 
♦P3*U+P4*Y2 
GOTO  1090 


1105  IF  I#0  THEN  1120  ELSE  1=1 

1110  MOVE  U+P1,S+P3  0  DRftW  U-Pl, 
S-P3  0  MOVE  U+P2.S-P4  ft  DRfl 
M  U-P2.S+P4 

1115  I8=SQR<1 . 39'L9>  0  T=T*I8  0 
Y 1 =Y 1216  0  GOTO  1075  0  11. 
39  IS  50%  CHI -SQUARED  VALUE 

1120  MOVE  K3- . 225TS 1 , K6+ . 1 TS2  0 
LABEL  "50%  &  95%  CONFIDENCE 
ELLIPSES" 

1125  X=LEN< "MU="&Bf > 

1130  MOVE  K3- . 925TXTS1 >  K5- . 22S2 
0  LftBEL  MMU= “&B$ 

1135  Y=LEN<  "SIGMft="lrft$> 

1149  MOVE  K4+< . 1-  05*<Y-1>>*S1,K 
5-.2TS2  0  LABEL  "SIGMA="fcA* 

1145  MOVE  U+ .  6*S 1 >  S- . 562S2  0  LAB 
EL  "U“  0  MOVE  K3,K6+.02*S2 
0  LABEL  "S"  0  IF  P8>=2  THEN 
COPY 

1150  RETURN 

1155  1  SUBROUTINE  OF  FOR  SCALING 
THE  N0 . ’ S  MU  AND  SIGMA  ON 
GRAPH 

1160  2=ABS<21>  0  IF  FPCZl )=0  AND 
ABS<Z1> <1000000080  THEN  AT 
=VAL*<21>  0  RETURN 

1165  IF  2<1  THEN  C*="-"  ELSE  C#= 

1170  T=LGT<2)  0  T 1 =FP<T> 

1175  IF  T 1 =8  THEN  1200  ELSE  18=1 
NT<T> 

1130  Tl=  5*<1-SGN<T>>+Tl  0  Tl=10 
-T1+. 000005 

1185  A*=VAL*<T1>  0  Z=SGN<21 >*VftL 
<A2C 1 / 73>  0  IF  FP<2>=0  THEN 
1205 

1190  A*=VAL*<Z>  0  IF  IO#0  THEN  ft 
*=A*8,C*&VftL*<ftBS<I8>> 

1195  RETURN 

1290  R*=VAL*<SGN<Z1>>8."0"  0  ft*= 
fl$8,C**-VftL$<IP<ABS<:T>>>  0  RE 
TURN 

1205  A*=VflL$<Z>S."  .0"«.CfiVftL$<ftBS 
< IP< I8>  >  >  0  RETURN 
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